{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "59b635f6-4d9f-4f95-bff8-4287d943f91a",
   "metadata": {},
   "source": [
    "# Elastic properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4726e73e-94b6-4d84-96dd-39253fdb0bbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ase.build import bulk\n",
    "from ase.calculators.emt import EMT\n",
    "from ase.io import read\n",
    "from ase.units import GPa, kg\n",
    "from calorine.calculators import CPUNEP\n",
    "from calorine.tools import get_elastic_stiffness_tensor, get_force_constants, relax_structure\n",
    "from matplotlib import pyplot as plt\n",
    "from pandas import DataFrame"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43ba3fde-7ecc-436a-b0b6-510bf9bdf333",
   "metadata": {},
   "source": [
    "This tutorial illustrates the calculation of the elastic stiffness tensor $c_{ij}$ and its relation to the sound velocity.\n",
    "In the first part we use silver in its face-centered cubic (FCC) ground state structure described by an effective medium theory (EMT) potential as a particular simple example.\n",
    "In the second part we then consider the orthorhombic structure of CsPbI<sub>3</sub>, which allows us to demonstrate the difference between relaxed $c_{ij}$ and clamped elastic constants $c_{ij}^0$.\n",
    "\n",
    "For background on crystal elasticity and the role of symmetry you can consult, e.g., Nye, *Physical Properties of Crystals: Their Representation by Tensors and Matrices*, Oxford University Press (1957)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4406f9a1",
   "metadata": {},
   "source": [
    "All models and structures required for running this and the other tutorial notebooks can be obtained from [Zenodo](https://zenodo.org/record/10658778).\n",
    "The files are also available in the `tutorials/` folder in the [GitLab repository](https://gitlab.com/materials-modeling/calorine/-/tree/master/tutorials)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa3f7c76-8e88-4b4f-86ea-547e67660845",
   "metadata": {},
   "source": [
    "## Face-centered cubic silver\n",
    "\n",
    "### Elastic constants and relation to sound velocities\n",
    "\n",
    "First we set up the primitive structure of FCC Ag and relax the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "95bb7ece-9e8f-4a9a-af69-03533ed14e16",
   "metadata": {},
   "outputs": [],
   "source": [
    "structure = bulk('Ag')\n",
    "calculator = EMT()\n",
    "structure.calc = calculator\n",
    "relax_structure(structure, fmax=0.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8fa677a-303c-452a-84ca-f943ea96b219",
   "metadata": {},
   "source": [
    "The elastic stiffness tensor $c_{ij}$ can then be readily calculated using the `get_elastic_stiffness_tensor` function.\n",
    "The latter applies a series of deformations to the cell and fits the resulting strain energy to a Taylor expansion in strain, in which the components of the elastic stiffness tensor appear as expansion coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "86b5bb00-2e5b-42ae-9f43-3203e5287c20",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[125.9  87.3  87.3   0.   -0.    0. ]\n",
      " [ 87.3 125.9  87.3  -0.    0.    0. ]\n",
      " [ 87.3  87.3 125.9  -0.   -0.    0. ]\n",
      " [ -0.   -0.   -0.   54.2  -0.   -0. ]\n",
      " [ -0.    0.   -0.   -0.   54.2  -0. ]\n",
      " [  0.    0.    0.   -0.   -0.   54.2]]\n"
     ]
    }
   ],
   "source": [
    "cij = get_elastic_stiffness_tensor(structure)\n",
    "with np.printoptions(precision=1, suppress=True):\n",
    "    print(cij)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "814feee6-0c3f-4faa-830c-8d946387754a",
   "metadata": {},
   "source": [
    "The elastic constants are related to the sound velocity along different crystal directions.\n",
    "In particular $c_{11}$ is related to the longitudinal speed of sound along $\\left<100\\right>$ in the long-wavelength limit (i.e., $\\boldsymbol{q}\\rightarrow 0$) according to\n",
    "\n",
    "$$\n",
    "c_s^{\\left<100\\right>} = \\sqrt{\\frac{c_{11}}{\\rho}},\n",
    "$$\n",
    "\n",
    "where $\\rho$ is the mass density of the material.\n",
    "\n",
    "Let us now test this relationship.\n",
    "For simplicity of presentation, we use the conventional 4-atom unit cell in the following, for which the cell axes are oriented along the Cartesian directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "23385494-8037-4615-a2fd-f0b409e386f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "structure = bulk('Ag', cubic=True)\n",
    "calculator = EMT()\n",
    "structure.calc = calculator\n",
    "relax_structure(structure, fmax=0.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07274b4b-73ec-4b5c-ace7-21ee662687bd",
   "metadata": {},
   "source": [
    "We compute the force constants using the `get_force_constants` function and then use methods of the [Phonopy object](https://phonopy.github.io/phonopy/phonopy-module.html) returned by the latter to calculate the group velocities at $\\boldsymbol{q}=(\\delta, 0, 0)$ with $\\delta = 0.01$.\n",
    "This $\\boldsymbol{q}$-point lies along $\\left<100\\right>$ and the small $\\delta$-value mimics the long-wavelength limit, i.e., $\\boldsymbol{q} \\rightarrow 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0102bb7b-e7b2-44be-8d85-a7b55eb349a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "phonon = get_force_constants(structure, calculator, [2, 2, 2])\n",
    "phonon.run_band_structure([[[0.01, 0, 0]]], with_group_velocities=True)\n",
    "band = phonon.get_band_structure_dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1f05bf4-07af-4b80-b12a-0f06d0b73cba",
   "metadata": {},
   "source": [
    "Next we select the speed of sound of the LA branch and convert it to SI unitss (m/s)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "195047a5-b2af-455c-b7b9-3d2a40165162",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "speed of sound, LA, 100 :  3415.1 m/s\n"
     ]
    }
   ],
   "source": [
    "group_velocities = np.linalg.norm(band['group_velocities'][0], axis=2) * 1e-10 / 1e-12  # Å/ps --> m/s\n",
    "speed_of_sound = group_velocities[0][2]\n",
    "print(f'speed of sound, LA, 100 : {speed_of_sound:7.1f} m/s')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29290fb1-c931-4dd5-9a39-5beb641860c7",
   "metadata": {},
   "source": [
    "Combining the speed of sound with the mass density, then yields the elastic constant $c_{11}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5205ea40-60bf-4681-8da1-96e7227fce06",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "density                 : 10677.9 kg/m^3\n",
      "elastic constant c11    :   124.5 GPa\n"
     ]
    }
   ],
   "source": [
    "density = np.sum(structure.get_masses()) / structure.get_volume() / kg / 1e-30  # amu/Å^3 --> kg/m^3\n",
    "print(f'density                 : {density:7.1f} kg/m^3')\n",
    "\n",
    "elastic_constant = speed_of_sound ** 2 * density * 1e-9  # Pa --> GPa\n",
    "print(f'elastic constant c11    : {elastic_constant:7.1f} GPa')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06dc0dcf-deb3-4123-b6c7-3986ee6d5de4",
   "metadata": {},
   "source": [
    "The result of 124.5 GPa is in good agreement with the value of 125.9 GPa obtained above by direct evaluation of the elastic stiffness tensor."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b5d212e-2bb8-4026-8166-cadb2ade8f03",
   "metadata": {},
   "source": [
    "### Presssure dependence\n",
    "\n",
    "We can now also readily evaluate the variation of the elastic constants with pressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c0d887a1-53b0-4887-8196-ed2094e797a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = []\n",
    "for volsc in np.arange(0.88, 1.05, 0.04):\n",
    "    structure_strained = structure.copy()\n",
    "    cell = structure_strained.get_cell()\n",
    "    cell *= volsc ** (1 / 3)\n",
    "    structure_strained.set_cell(cell, scale_atoms=True)\n",
    "    structure_strained.calc = calculator\n",
    "    relax_structure(structure_strained, constant_volume=True)\n",
    "    \n",
    "    cij_rlx = get_elastic_stiffness_tensor(structure_strained)\n",
    "    pressure = -np.sum(structure_strained.get_stress()[:3]) / 3 / GPa\n",
    "    \n",
    "    data.append(dict(volsc=volsc,\n",
    "                     pressure=pressure,\n",
    "                     volume=structure_strained.get_volume() / len(structure_strained),\n",
    "                     c11=cij_rlx[0][0],\n",
    "                     c12=cij_rlx[0][1],\n",
    "                     c44=cij_rlx[3][3],\n",
    "                    ))\n",
    "df = DataFrame(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ae9d8c4-62fa-4ff1-a2e5-68d21c148aee",
   "metadata": {},
   "source": [
    "The results show a stiffening of the elastic constants with pressure, which is expected for this system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "664bc3fd-c940-4d46-b237-ae3c4f4bb4ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFeCAYAAABKNlxZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAABWIAAAViAHE10CgAAB2mElEQVR4nO3deVhUZfsH8O8AAwwMMOwMsivK6oKAosjinlZKqS0/d00tzbe0PRcy7bXXLH3L0jLTbDHK9C3L3GIRxAVZREgFWWRHtmEdGGbO74+RkXFYBhiYAe7PdXEV5zznnOcQzc2z3Q+LYRgGhBBCCNFIWuquACGEEELaR4GaEEII0WAUqAkhhBANRoGaEEII0WAUqAkhhBANRoGaEEII0WAUqAkhhBANRoGaEEII0WAUqAkhhBANRoGaEEII0WAUqAkhhBANRoGaEEII0WA66q7AYOLm5obKykq4uLiouyqEEEL6UFZWFkxNTXHr1q0uX0uBug9VVlaivr5e3dUghBDSx3ry2U+Bug+1tKTj4+PVXBNCCCF9KSAgoNvX0hg1IYQQosEoUBNCCCEajAI1IYQQosFojJoQQghph1Akxo18AWobReDqsTHSzgT6bO0+rQMFakIIIeQRxQIhvo7NQkRCPgQNItlxEw4bC3ztsHKSC6yN9fukLhSoCSGEkFbSCgVYcugqymqbFM4JGkT46mI2TiQV4Mhyf3jamvR6fWiMmhBCCHmgWCBsN0i3VlbbhCWHrqKkWtjrdaJATQghhDzwdWxWp0G6RVltE76Oze7lGlGgJoQQQgBIJ45FJOR36ZqIa3kQisS9VCMpCtSEEEIIgBv5ArmJY8qoahAhtUDQSzWSokBNCCFk0GtqliD6dmm3rq0Rdi24dxXN+iaEEDJo3S6uQURCHk4kFaCiTrmx6UcZ6bNVXCt5FKgJIYQMKoIGEX5PKcTPCXlIye9ZtzWPw4b3kN5dokWBmhBCyIAnkTC4nFWOiIQ8nL5ZjMZmiUIZryHGMNZn49LdcqXvu8DPvtczlVGgJoQQMmAVVDXgl4R8/Hw9D/mVDQrneQZszB09BPN97eBpa4KSaiFm//eiUku0LLl6WBHo3BvVlkOBmhBCyIAiFIlxNr0EPyfkITazDAwjf57FAia5WuIZX3tM9bCCns7DFrG1sT6OLPfvNOmJJVcPR5b790kaUQrUhBBC+j2GYZBWWI2IhDycTCpAtbBZoYyDmQHmj7XD02PtYMvjtHsvT1sT/LF+Er6OzUbEtTxUtVqyxeOwscDPHisCnSnXN+ldmrAjTE/V19fjiy++wPXr15GYmIg7d+6AYRhkZ2fDycmpzWt+//13nD17FomJiUhOTkZ9fT22bt2K8PDwPq07IUQ1KuuacDK5ABEJ+finqFrhvD5bC7O8+Vjgaw9/JzNoabGUuq+1sT7emeWODdOGI7VAgBqhCEb60oljtHsW6VWatCNMT5WWluK1114DADg6OsLU1BQVFRUdXrN7925ER0fD2NgYtra2yMzM7IuqEkJUSCxhcDHjPn5OyMe59BI0iRUnho225+EZP3vMHsmHcQ+WT+mzteHnZNaT6vYYBepBRNN2hOkpCwsLnD17FmPHjoWZmRlmzpyJM2fOdHjN+++/DxsbGwwbNgw//fQTnnvuuT6qLSGkp3LL6/BzQj6OJ+ajSKC4GYa5oS6e8hmC+b72GG5tpIYa9g4K1INEV3eE+WP9JLW3rGNiYvDxxx8jPj4eVVVVsLKygp+fHzZs2IDAwEBwuVxMmzatS/ecNGlSL9WWENIb6puacTq1GBEJebiSrdhjpq3FQugIKyzwtUOomxXY2gMv4SYF6kGiOzvCvDPLvZdr1b69e/fi1VdfBYfDQVhYGBwcHFBQUIDY2Fj88ssvCAwMVFvdCCG9i2EYJOVV4eeEPPyeUoTaRsWJYS6Whljga4+nfIbAyqh/DNd1FwXqQaC7O8JsmDZcLRPMUlJSsGHDBvD5fMTFxclNDGMYBkVFRX1eJ0JI77tf04gTSfmISMhHZmmtwnlDXW08PtIWC/zs4ONgChZLuYlh/V2/CdQFBQX4+eef8eeff+LWrVsoLi6GmZkZJk6ciDfeeAPjxo1TuKa6uhrh4eE4fvw4iouLwefzMX/+fGzduhVcLlehvEQiwb59+/Dll18iMzMTXC4XU6dOxY4dO+Di4tIXr9kpn/fPoaGpa1uqiSVMm5MtOlLVIMLI8LPQVnKGZGscXW0kbu5al3RrBw4cgEQiwfbt2xVmb7NYLNja2nb73oQQzdIsliDy9n1EJOQh8lYpmiWMQhl/JzPM97XDLG8+DPX6TdhSmX7zxp9++ik+/PBDDB06FNOnT4elpSUyMjJw8uRJnDx5Ej/88AOeeeYZWfm6ujoEBwcjOTkZ06dPx3PPPYekpCR89NFHiI6ORkxMDPT15btLVq9ejYMHD8LT0xPr169HYWEhIiIicPbsWVy+fBmurq59/doKGprEaOjlvU9bNIklQN88Ss7Vq1cBANOnT+/7hxNC+kRmaS1+TsjDr0kFuF/TqHDe2lgPT/vYYd5YO7hYKjasBhOVBeqsrCz8/fffiIuLQ35+PsrKymBgYABLS0t4e3sjODgYQUFB0NXV7db9/f39ERUVheDgYLnjFy9exJQpU/Diiy9i7ty50NPTAwD85z//QXJyMt58803s3LlTVv6tt97Chx9+iE8++QRvv/227HhkZCQOHjyIoKAgnDt3TlbP559/HrNmzcK6des6nVHcFzi6Xe+K7k6LGgB0tbW63aLuCYFAABaLBT6f36P7EEI0S41QhD9uFCEiIQ+J96oUzrO1WZjqbo0FvvaY5GoBnQE4Maw7ehSoGYbBsWPHsH//fsTGxsqOPeq3337DBx98AFNTUyxduhRr166Fs3PX8qM+9dRTbR6fNGkSQkNDcfbsWaSmpsLX1xcMw+DgwYPgcrnYvHmzXPnNmzdj3759OHjwoFyg/uqrrwBIl++0/mPiscceQ0hICM6ePYt79+7BwcGhS/VWte50KQtFYoz74EKXNkTncdi4/M4UtYxR83g82Vj0kCFD+vz5hBDVYRgG13Iq8dO1PPyZWtRmj+AIayMs8LPH3NG2MOfqqaGWmq3bf6789ddfGDVqFP7v//4P//zzD1asWIGDBw8iJSUFxcXFaGpqgkAgQHZ2Nv766y+Eh4fD3d0dn3zyCdzd3bFhwwZUVlaq5CXYbOlidh0d6d8dGRkZKCwsxMSJE2FoaChX1tDQEBMnTkRWVhby8vJkx6OiomTnHjVjxgwAQHR0tErq29f02dpY4GvXpWv6YkeY9vj7+wMAzp49q5bnE0J6rlggxL7ITIR+FIUFB+JxPDFfLkgb6etg4XgH/LZuIv56ZRJWBDpTkG5Ht1vUs2bNQmBgIH777TfMnDlTFiRbMzIygpGRERwdHTF9+nRs3rwZubm5+Oqrr/DZZ5+Bx+Nhy5YtPXqBe/fu4fz58+Dz+fD29gYgDdQA2h1TdnV1xZkzZ5CRkQF7e3vU1dWhqKgIXl5e0NZWDE4t92m5b3+0cpILTiQVaNSOMO1Zs2YNDhw4gE2bNmHy5MlwdHSUnWtpadOEMkI0T1OzBBf+KUFEQh6i79xHG/PCMHGYORb42mOGp02/S1usLt0O1OfOncOUKVO6fJ2joyO2b9+O1157DdnZ2d19PABAJBJh0aJFaGxsxIcffigLsgKBdCNwE5O2s2sZGxvLletq+c4EBAS0efzmzZvw8vJS6h6qpok7wrTH29sbe/bswfr16+Hp6Ym5c+fC0dERxcXFiImJwezZs7Fnzx4AwGuvvYaysjIAQGpqquxYy6z+lStXyq25bpl8CED2+3fy5Enk5OQAANzc3PDWW2/1wVsSMnDcKq5GxLV8nEwuQEWd4ufLEB4H88ZKJ4bZmxmooYb9W7cDdXeCdGs8Hg9jxozp9vUSiQRLly5FTEwMXnjhBSxatKhH9RkMNG1HmI6sW7cOXl5e2L17N06fPo3a2lpYWVlh3LhxWLBggazcL7/8gtzcXLlrjx8/Lvv3kJAQuUCdnJyMI0eOyJVPSUlBSkoKACA4OJgCNSFKEDSI8FtKIX5OyMONfMVGjK6OFmZ62mCBrz0mDDVXejMMoqjfLM9qTSKRYPny5fjhhx+wcOFC7N+/X+58S8u4vRZwdXW1XLmulu9MfHx8m8fba2n3JU3aEaYzISEhCAkJ6bBMS0tYWeHh4bRTFiHdJJEwiM8qR0RCHv66WYzGZsXVJN5DTLDA1w5PjhoCE4Pub4ZBHuqVQC0Wi1FWVobGRsW1cQB6NHNaIpFg2bJl+Pbbb/Hcc8/h8OHD0NKSnxPX2Zjyo2PYhoaG4PP5yM7OhlgsVhin7mzMuz/ShB1hCCH9Q35lPX65no+fE/JRUNWgcN7UgI25Y4Zg/lh7eNgaq6GGA5tKA/X169fxzjvvICYmBk1NbY+DslgsNDcr5m1VRusg/cwzz+Do0aPtTv6ytbVFXFwc6urq5GZ+19XVIS4uDs7OzrC3t5cdDw4OxrFjxxAXF4egoCC5+7Wsn370OCGEDFRCkRhn0orxc0I+4u6W4dGVt1osIGi4JRb42mOKuxX0dDSrR24gUdlq8uTkZEyaNAnx8fGYPn06GIbByJEjMX36dFhYWIBhGAQHB3d7LLmlu/vbb7/F/Pnz8d1337UZpAHpHwMrV65EbW0t3n//fblz77//Pmpra/HCCy/IHV+1ahUA6Trr1n9knD59GlFRUZg+fbrc7GNCCBloGIZBar4Am0/ehP+O8/jXsWTEZsoHaUdzA7w+YwTi3pqMw8v8McubT0G6l7GYtjKUdMPTTz+N06dP4/r163B3d4eWlhbCw8OxZcsWNDQ0YOPGjfjll19w9epVhfzNyggPD8d7770HLpeLf/3rX20uB5s7dy5Gjx4NQNpynjhxIlJSUjB9+nT4+PggMTERZ8+ehZ+fH6Kjo8HhcOSuf+GFF2QpRGfPno2ioiL89NNP4HK5iI+Px/Dhw7vzo5FpGaNubwybEELUoaKuCSeTChCRkIdbxTUK5zlsbczy5mOBrx38nc0GzWYYqtSTz3+VdX3HxsbiySefhLv7w60RW/4G4HA4+Oyzz3Dp0iW88847+OGHH7p8/5ZJQ7W1tdixY0ebZZycnGSB2tDQENHR0bJNOSIjI8Hn87Fx40Zs3bpVIUgD0s0gvL298eWXX2Lv3r3gcrkICwvDjh07MHTo0C7XmRBCNJVYwiAm4z5+TsjDufQSiMSKbTYfBx4W+Npj9kg+jPRpYpi6qCxQCwQCuR2m2Gw2amsfblOmpaWFkJAQ/Pjjj926/+HDh3H48OEuXWNiYoJPPvkEn3zyiVLltbS0sH79eqxfv74bNSSEEM2XU1aHn6/n4fj1AhRXCxXOW3B18bSPHeb72mGYlZEaakgepbJAbWVlJZcS1MbGRmHWtVAoRH19vaoeSQghRAn1Tc34M7UYEQl5uJpdoXBeW4uFyW5WWOBrj5ARlmDTZhgaRWWB2sPDA7dv35Z9P3HiRJw8eRLx8fEICAjAP//8g4iICLi5uanqkYQQQtrBMAwS71Xh54Q8/J5SiLo29rEfammIZ/zsMXfMEFgZqT/REWmbygL17Nmz8eqrr6KoqAh8Ph9vvvkmTpw4gcDAQJiZmaGyshISiQTvvPOOqh5JCCHkEaU1QpxIlE4Mu3u/TuE8V08HT4ziY95Ye/g48GhiWD+gslnfIpEIFRUVMDU1lW0TeenSJezYsQNZWVlwdHTEyy+/jNmzZ6vicf0SzfomhPQGkViCyFuliEjIR+TtUojb2A3D39kMC3ztMcvbBga6/TIpZb+m9lnfQqEQaWlpYLFY4PF4suMTJkzAH3/8oYpHEEIIeURmaQ0iEvLxa2J+m5vtWBvrPdgMwx7OFoZt3IH0Bz0O1J988gm2bNkimyRmaGiI7du308xpQgjpBTVCEU7dKEJEQh6S7lUpnGdrszDNwxrzfe0R5GoJbdoMo9/rUaD+7bffsHHjRrBYLIwYMQIAcPv2bbz66qsYOnTooO7mJoQQVWEYBleyKxCRkIc/U4sgFCluhuFmY4QFvtKJYWaGumqoJektPQrU+/btg46ODv78809MnToVAHDhwgU89thj2LdvHwVqTSZqAAqTgMYaQM8IsB0DsBWTwGiy+vp6fPHFF7h+/ToSExNx584dMAyD7OzsNrPflZeX4/jx4zh16hRu3ryJgoICGBkZwc/PD6+88gpmzJjR9y9BSAeKBA04fj0fP1/PR2654tJWY30dzBk9BAt87eE1xJgmhg1QPQrUiYmJePLJJ2VBGpDuUz1nzhxER0f3uHKkF1QXAvH7gKTvAGHVw+P6PGDMQiBgHWDMV1ftuqS0tBSvvfYaAMDR0RGmpqaoqFBcI9ri559/xosvvghbW1tMmTIFQ4YMQX5+Po4fP46//voL//nPf/D666/3VfUJaVNjsxjn00sRkZCHixn30ca8MAQOs8B8XzvM8LTRuO1pier1KFBXVla2uS56xIgROHnyZE9uTXpD0Q3gu6eBulLFc8IqIP4z4EYEsPA4wB/Z59XrKgsLC5w9exZjx46FmZkZZs6cKdvprC3Dhw/Hb7/9htmzZ8ttjbpp0yaMGzcO7777Lv7v//4Ptra2fVF9QuSkF1YjIiEPJ5MLUFUvUjg/hMfBfF87PO1jB3szAzXUkKhLjwK1RCKRLcVqjc1mQyJRHEMhalRd2H6Qbq2uVFpudYzaW9YxMTH4+OOPER8fj6qqKlhZWcHPzw8bNmxAYGAguFwupk2bpvT9Jk+e3ObxESNG4JlnnsGXX36JS5cuYd68eap6BTJICEVi3MgXoLZRBK4eGyPtTJRq6QrqRfhfinTN882CaoXzujpaeMzLBgt87RHgYg4tmhg2KPV41jeNifQT8fs6D9It6kqBy/uA6dt7t04d2Lt3L1599VVwOByEhYXBwcEBBQUFiI2NxS+//ILAwECVPo/Nlm440NaubIS0p1ggxNexWYhIyIeg4WEr2ITDxgJfO6yc5AJrY/mMXxIJg7i7ZYhIyMeZtGI0NSs2akbamWC+rz2eHGkLEwPaDGOw61HCEy0tLejo6Ch8uDU3N0MsFkNPT0/xgSwW6uoUs+UMBmpLeCJqAHa7yY9Jd4ZjCmz4Ry0TzFJSUuDj4wMbGxvExcXJTQxjGAZFRUVtdk+3dH23N5msPdXV1Rg+fDgEAgHy8/Nhbm6ugrcgA11aoQBLDl1tc/1yCwuuLo4s94enrQnyKurx8/V8HL+ej4KqBoWypgZshI2RbobhzjfuzaoTNVBbwhMHBwdqUfe1/7hIA29XSMSAuLFr1zRUAjsdAa1uTFRhc4A3srp+3QMHDhyARCLB9u3bFQIui8VS+RjymjVrUFJSgm3btlGQJkopFgg7DdIAUFbbhGcPXMYIGyMk5FYqnNdiAcHDLbHA1x5T3K2hq0ObYRBFPQrULXtEkz4kagBEfbQDmbgRUMzj3+uuXr0KAJg+fXqvP+vtt9/Gjz/+iJkzZ1IeeqK0r2OzOg3SLWoamxWCtJO5Aeb72uNpHzvYmNBmGKRjNCDX33SnK7o7LWoA0Nbrfou6BwQCAVgsFvj83p3MtnnzZuzcuROTJ0/Gr7/+Cm1tWuZCOicUiRGRkN/l6/R1tPD4KFss8LWHn5Mp9UYSpVGg7m+606Xcz8aoeTyebCx6yJAhvfKMzZs3Y/v27QgJCcHvv/8ODqd/JXsh6nMjXyA3cUxZXy72RdBwy16oERnoehyor1+/jurqagQHB8vWpv7vf//D//73P4Wy48aNw+rVq3v6SNJVbI40mUn8Z8pfM2ah2jKV+fv7IyEhAWfPnsWyZctUfv+WIB0cHIw//vgDBga0JpUor7ax60EaAJppySrpph4F6nv37mHixIl46qmnEBoaKjuenJyMw4cPK5T/4YcfMGvWLNjb2/fksaQ7AtZJk5kos0SLaw2MX9v7dWrHmjVrcODAAWzatAmTJ0+Go6Oj7FxHs76VsWXLFmzfvh2TJk2iIE26rK6xGbGZZd261kifllmR7ulRoD5y5Aiam5uxbds2hXMsFguXLl2SfV9ZWYnZs2fj8OHD2Lx5c08eS7rDmC/NONZZ0hOutbScGpOdeHt7Y8+ePVi/fj08PT0xd+5cODo6ori4GDExMZg9ezb27NkDAHjttddQVib94ExNTZUd43K5AICVK1fK1lwfPnwY77//PnR0dODv749du3YpPDskJAQhISG9/5KkXymoasCRSzn48eo91Aibu3w9j8OG9xCTXqgZGQx6FKjPnz+PsWPHYtiwYW2eHzdunNz3QUFBOH/+PAVqdeGPlGYcu/wg13dDq5moHFNpd/f4tWrPSAYA69atg5eXF3bv3o3Tp0+jtrYWVlZWGDduHBYsWCAr98svvyA3N1fu2uPHj8v+PSQkRBaoW1YpNDc3Y/fu3e0+mwI1AaS9N4n3qnAoNht/pRVD3FbSbSUt8LOnnNyk23oUqP/55x/Mnz9f6fJubm5yH6JEDYz50oxjoe8ChclAYzWgZwzYjta43bOUad12ZYlgeHg4wsPDe1QnMvCJxBL8mVqEQ3E5SMmrUjjv62iKp3yG4ONzd5RaomXJ1cOKQOdeqCkZLHoUqKuqqmBmZqZwfO7cuW1mhjI3N4dAIOjJI4mqsDmAY4C6a0GIxqiqb8IPV+/h20u5KK4Wyp3T0WJh9kg+lk10xmh7HgBglD2v06Qnllw9HFnur5BGlJCu6FGg5nK5qK5WTCQ/atQojBo1SuF4dXU1Td4hhGiUzNJafBOXjeOJ+RCK5Gdm8wzYeN7fAYsCHME3ke9x8rQ1wR/rJ+Hr2GxEXMtDVaslWzwOGwv87LEi0JmCNOmxHgVqR0dHJCQkKF3++vXrXcrBTAghvYFhGFzMKMOhuGxE3b6vcH6opSGWBzrjqTF24Oi2P7ZsbayPd2a5Y8O04UgtEKBGKIKRvnTiGI1JE1XpUaAODQ3F3r17kZSUhDFjxnRYNjExEfHx8Xj11Vd78khCCOk2oUiME0kFOBSbjYzSWoXzQcMtsXyiE4JcLbu0paQ+Wxt+TorDgISoQo8C9bp16/DZZ59hwYIFOH36dLuzv+/evYtnnnkGOjo6WLtWfetzCSGDU0m1EEfjc/H9lVxU1ssnLNHT0cJTPnZYPtEJrtZGaqohIe3rUaB2cXHBzp078dprr2HkyJFYsGABQkJCZMkoCgsLERkZiZ9//hlCoRC7d++Gi4uLSipOCCGdSc0X4FBcNk7dKIRILL+8yspID0smOOE5fweYGeqqqYaEdK7HKUQ3bNgAQ0NDvPnmm/j2229x9OhRufMMw8DIyAh79uzBqlWrevo4QgjpkFjC4Fx6MQ7F5uBqToXCee8hJlgR6IxZ3nzaVpL0CyrZlGP16tV45pln8MsvvyA+Ph7FxcUAAGtra0yYMAHz5s0Dj8dTxaMIIaRN1UIRIq7l4fClHORXyu/ZrsUCZnjaYEWgM8Y60s5VpH9R2e5ZPB4PK1euxMqVK1V1S0II6VRueR0OX8rBzwn5qG2UT+9ppKeDZ/zssWSCE+zNaGko6Z9om0tCSL/DMAyuZFfgUGw2zv1TAuaR7J6O5gZYNsEJ83ztwdWjjznSv3X7N/jYsWN49tlnu/3gvLw82e5bhBCijMZmMU6lFOFQXDbSChWTLY13McOKQBdMdrOCdheWVxGiybo9k2LhwoUYOXIkjhw5gtpaxfWI7UlISMDq1avh6uqKyMjI7j6eEDKIlNU24r8XMhD4YSQ2/pwiF6R1tbXwtI8d/lgfiGOrAjDNw5qCNBlQut2ivnbtGjZs2IBly5Zh7dq1mDlzJsaNG4exY8fC2toaPB4PQqEQFRUVuH37Nq5cuYJz584hIyMDxsbG2LZtG/71r3+p8l0IIQPMreJqHIrNxsnkQjQ1y6f3NDfUxcLxjvi/8Q6wMqI0nWTg6nagHjNmDCIjI/H333/jiy++wG+//YZff/213dmUDMPI1l2vXLkSpqam3a40IWTgkkgYRN0pxdex2YjLLFc472ZjhOWBznhylC2l6SSDQo9nWUyePBmTJ09GVVUVYmJiEBcXh/z8fJSXl4PD4cDS0hLe3t4IDg7GyJEjVVFnogLCZiHSytNQJ6qDIdsQnuae0Nfp/62SWbNm4fTp09DT04NQKOy0/E8//SSba/Hjjz/2aN4F6Zm6xmYcT8zHN3E5yC6rkzvHYgFT3KywfKIzAoaa0/IqMqiodHnWk08+iSeffFJVtyS9oKSuBEfTj+JE5glUNz0c5zPWNUbYsDAs9lwMKwMrNdaw+7766iucOXMG+vr6YB6dBtyG4uJirF27FoaGhqirq+u0POkdBVUN+PZSDn68eg/VQvnlVQa62pg/1g5LJzrD2cJQTTUkRL1o3cIgcqviFtacW4NyoWJ3YnVTNY6kH8GprFPYP20/3Mzc1FDD7svJycHGjRuxYcMG/Pzzz7KkOx1ZtWoVjIyMsHTpUuzevbsPaklaS7xXia9js/HXzWKIJfJ/WA3hcbBkgiOe8XWAiQFbTTUkRDNQ/rxBoqSupN0g3Vq5sBxrzq1BaX1pH9WsfTExMZg7dy6sra2hp6cHe3t7PPXUU4iNjZUrxzAMli9fDj6fj23btil178OHD+P333/HwYMHweVye6P6pA0isQS/pRRi7r44PPX5Jfxxo0guSI91NMW+530Q/XoIVgUNpSBNCKhFPWgcTT/aaZBuUS4sx9H0o9jou7GXa9W+vXv34tVXXwWHw0FYWBgcHBxQUFCA2NhY/PLLLwgMDJSV/fTTTxEdHY2YmBhwOJxO752Xl4dXXnkFq1atwpQpU3Dx4sXefBUCoKq+CT9ezcO38TkoEsjPHdDRYmGWNx/LA50x2p6nngoSosEoUA8CwmYhTmSe6NI1JzJOYO3otWqZYJaSkoINGzaAz+cjLi4OTk5OsnMMw6CoqEj2fUZGBt5++22sX79eqeQ5DMNgxYoVMDY2xkcffdQb1Set3L1fi2/isnH8egEaRGK5cyYcNp4f54DFAY7gm3T+BxYhgxUF6n4m6FgQhOLOZzO3JmbEaBI3dekaQZMAE49NhDar68tf9LX1EfNsTJeva3HgwAFIJBJs375dLkgDAIvFkm2jKpFIsGTJEvD5fOzYsUOpe+/fvx/nzp3DX3/9BSMj2nu4NzAMg9jMMhyKzUbk7fsK510sDbF8ojOe8hkCA136CCKkM/R/ST8jFAvR0NzQeUEV6GpwV5WrV68CAKZPn95huV27duHy5cuIjIyEgUHnGy5kZWXh9ddfx/LlyzFjxgyV1JU8JBSJcTKpAIfisnGnRDFb4SRXC6wIdEaQqyW0KHMYIUqjQN3P6Gt3vSu6Oy1qANDV1u12i7onBAIBWCwW+Hx+u2Xu3LmDrVu34qWXXkJwcLBS912xYgV4PB4+/vjjHtWPyCutFuLo5Vx8f+UeKurkf8/0dLTwlI8dlk10wnBr6sEgpDtYjDILTpXg4uKCV155BevXr2+3zL59+7B7925kZWWp4pH9TkBAAAAgPj6+T58rbBZiys9T5NZNd8ZE1wTn559Xyxi1n58fEhISkJ+fjyFDhrRZ5uTJkwgLC1PqfpWVleDxeODxeBAIBJ2W/+STT/DKK690pcqD0s0CAQ7FZuP3G4UQieU/RqyM9LBkghOe83eAmaGummpIiOboyee/ylrUOTk5qKqq6rBMVVUVcnNzVfVIoiR9HX2EDQvDkfQjSl8T5hqmtkxl/v7+SEhIwNmzZ7Fs2bI2yzg5OWHFihVtnvvpp5/Q0NCApUuXAgD09PQAAIsXL0Z9fb1C+cTERCQlJSE0NBQuLi7w8vJSzYsMQGIJg3PpJTgUl42r2RUK572HmGBFoDNmefOhq0OrPwlRhT7t+hYIBLIPTdK3FnsuxqmsU0ot0bLgWGCRx6I+qFXb1qxZgwMHDmDTpk2YPHkyHB0dZedaZn2PHj0aBw8ebPP68+fPo7i4WOH8f//73zbLh4eHIykpCatWraIUou2oEYoQkZCPw5eykVchP0dCiwXM8LTB8kBn+DqaUnpPQlSsR4E6JkZ+Zm9OTo7CMQAQi8XIy8vD999/j+HDh/fkkaSbrAyssH/a/k6TnlhwLLB/6n61phH19vbGnj17sH79enh6emLu3LlwdHREcXExYmJiMHv2bOzZs0dt9RtM7pXX4/ClHEQk5KG2UT69p5GeDp7xs8eSCU6wN+t8Mh8hpHt6FKhDQkJkfz2zWCwcOXIER4603b3KMAxYLBZ27tzZk0eSHnAzc0PEExHSXN8ZJyBoejhea6JrgjDXMCzyWKQRub7XrVsHLy8v7N69G6dPn0ZtbS2srKwwbtw4LFiwQN3VG9AYhsHV7Ap8HZuNc/+U4NFZLA5mBlg20Qnzfe3B1aP5qIT0th5NJgsPDweLxQLDMNi2bRuCg4MREhKiUE5bWxtmZmYIDQ2Fu7t7T+rbr6lrMllbhM1CpJeno1ZUCy6bCw9zjwGxexbpvqZmCU7dKMTXsdlIK1SceDjexQzLJzpjirs1tGl5FSFdorbJZOHh4bJ/j46OxrJly7B48eKe3JL0EX0dffhY+6i7GkQDlNc24vsr93D0ci7u1zTKndPV1sITo2yxbKITvIaYqKmGhAxuKuu3ioyMVNWtCCHdIBSJcSNfgNpGEbh6bIy0M4E+u/118LeLa3AoNhsnkgvQ1CyRO2duqIv/G++IheMdYGVEPS2EqFOvDDDV1dWhqqoKYrG4zfMODg698VhCBqVigRBfx2YhIiEfggaR7LgJh40FvnZYOckF1sbSYCuRMIi+cx9fx2YjNrNM4V5uNkZYPtEZT4627TDIE0L6jkoD9ddff43du3fj9u3b7ZZhsVhobm5u9zwhRHlphQIsOXQVZbWKmecEDSJ8dTEbJ5IKcGDRWKQXVuObuBxkldUplJ3iZoUVgc4IGGpOy6sI0TAqC9RffPEF1q5dCx0dHQQFBcHOzg46OqptsH/33Xe4ePEirl+/jtTUVDQ1NeGbb76RJbZoLTw8HO+9916798rOzlbY8AEAzpw5gw8++ACJiYlgsVgYO3YsNm3ahClTpqjwTQjpuWKBsN0g3VpZbRPmfRGPR2eNctjamO9rh6UTnOBiSXtyE6KpVBZJ9+zZAwsLC8TGxvbaWulNmzYhNzcXFhYW4PP5SmU5W7JkSZsBmcfjKRz77rvvsGjRIlhaWsqC/08//YRp06YhIiIC8+bN6+EbEKI6X8dmdRqkW7QO0rYm+lgywQnP+jnAxIDdO5UjhKiMygJ1bm4uVq5c2asJTQ4ePAhXV1c4Ojpi586dePvttzu9ZunSpW0uGXtUZWUlXn75ZVhYWCAxMRF2dnYAgDfffBNjxozBiy++iBkzZtDWiEQjCEViRCTkd+kabS0WPpo3Ck+M4kNHm9J7EtJfqOz/Vj6f3+7kMVWZOnWqXDpJVfr5559RVVWFl19+WRakAcDOzg7r1q1DWVkZTpw40SvPJqSrbuQL5CaOKUMsYWBnxqEgTUg/o7L/Y5csWYLTp0+jrk5xooo6xcTE4MMPP8SuXbtw8uRJ1NYq7pMLAFFRUQDa3gO5Ze/i6OjoXqsnIV1R29i1IN2iRti96wgh6qOyru9NmzYhLS0N06ZNw86dO+Hj4wMuV/0TVLZu3Sr3PY/Hw969exUSs2RkZAAAXF1dFe7RcqylTGdaMtA86ubNm7QzE+mxirom/HGjqFvXGunTmDQh/Y3KAnXLrlgMwyA0NLTdcn21PGvUqFE4dOgQQkJCwOfzUVxcjFOnTmHLli1YunQpeDwennzySVn5ln2KTUwUsy8ZGxvLlSFEHQqrGvDVxSwcu5qHBlHXh5l4HDa8KbsYIf2OygL1pEmTNGr9ZVhYmNz3Tk5OWLduHdzd3TFt2jRs2rRJLlCrUnu5XNtraRPSkbv3a7E/6i5OJhdAJO52an4s8LOnJCaE9EMqC9QtY7yabsqUKRg6dChSU1NRXV0tay23tKQFAgHMzc3lrqmurpYrQ0hfSM0X4POoTPyVVqywg1XQcEs862ePLf+7qdQSLUuuHlYEOvdSTQkhvWlQ7lFnYWGBzMxM1NfXywK1q6srEhISkJGRoRCoOxq/JkSVGIZB/N1yfB51VyHFJ4sFzPLi48WQobINMhzNDTpNemLJ1cOR5f6yNKKEkP5l0AXquro6pKWlwdDQEBYWFrLjwcHB+PHHH3H27FmMHz9e7pozZ87IyhDSGyQSBuf+KcHnUXeRklcld46tzcJTY+ywOthFIYOYp60J/lg/CV/HZiPiWh6qWi3Z4nHYWOBnjxWBzhSkCenHVBqoxWIxIiIicP78eRQWFqKxsVGhDIvFwoULF1T5WAU1NTUoKipSSL7S0NCAF154ATU1NVi2bJlcitMFCxbgzTffxKefforly5fL1lLn5+fjs88+g4WFhcK4NyE9JRJL8FtyIfZH30VGqfzSQQNdbTzn74CVk5zBN+G0ew9rY328M8sdG6YNR2qBADVCEYz0pRPHaEyakP5PZYG6rq4O06dPx+XLl8EwDFgsFphWA2st3/dkwtnBgwcRGxsLAEhNTZUdaxkfDwwMxMqVK1FeXg43Nzf4+fnB3d0dNjY2KCkpwfnz55Gfnw9vb2/s2rVL7t6mpqb47LPPsGjRIvj4+OCZZ54BIE0hWl5ejp9++omykhGVaWgS46dr9/DVxWwUVDXIneMZsLF0ghOWBDjB1FBX6Xvqs7Xh52Sm6qoSQtRMZYF6+/btiI+Px7Zt2/DSSy/BwsIC4eHhWL16NWJiYvDOO+/Ax8cH33//fbefERsbiyNHjsgdi4uLQ1xcnOz7lStXwszMDC+99BKuXr2KP//8E5WVleBwOHB3d8f69euxbt06cDiKLZSFCxfCwsICH3zwAb755hu5TTmmTp3a7XoT0kLQIMLR+Bx8E5eD8jr5cWUbY32snOSM5/wdYKg36EalCCHtYDHMo/NJu2fEiBEwNzfHpUuXAABaWloIDw/Hli1bAEi7kEeNGoXXXntNqRzdA1HL8qz2lm+Rgau0RoivY7Px/eV7qG2UzyPgbGGINcEumDtmCPR0qKuakIGoJ5//Kvuz/d69e5g9e7bsey0tLbkxajs7O8yePRtHjhwZtIGaDD73yutxIOYufr6ej6Zmidw5T1tjvBQyDDO9bKCtpTk5CAghmkVlgdrQ0BBaWg9Th5uYmKCoSD7NoY2NDe7du6eqRxKisW4VV+OLqLv4PaUQkkf6rMY5m+Gl0GEIcrXQqCRBhBDNpLJA7ejoKBeEvby88Pfff6OxsRF6enpgGAYXLlwAn89X1SMJ0TgJORX4IuouLtwqVTg31d0KL4YMw1hHUzXUjBDSX6ksUE+ZMgXffPMNmpuboaOjgyVLlmDlypUICAjAlClTcOnSJSQnJ2Pjxo2qeiQhGoFhGETduY8vIu/iak6F3DltLRaeHGWLNcFDMcKGVg0QQrpOZYH6hRdegLm5Oe7fvw8+n4/ly5cjKSkJn3/+OZKTkwEATz/9NMLDw1X1SELUSixh8GdqEb6Iuov0omq5c7o6WnjG1x6rglxgb2agphoSQgYClc36bs/9+/eRlZUFR0dH2NjY9OajNB7N+h4YGpvF+DWxAAei7yKnvF7unJGeDhYGOGL5RGdYGumpqYaEEE2jMbO+eTyeLHd2C0tLS1haWgKQZgyrrKyEg4ODqh5LSJ+pbWzGj1fu4WBsFkqq5bPuWXB1sWyiMxYFOMKY9nwmhKiQVudFlOPs7Iy9e/d2WOa///0vnJ1pBx/Sv1TUNeHjc3cwceff2PHnP3JBegiPg21zPBH75mSsDR1GQZoQonIqa1EzDIPOetF7uZedEJUqrGrAVxezcOxqHhpEYrlzrlZcvBgyFE+MsgVbW2V/7xJCNI2oAShMAhprAD0jwHYMwG4/935v6NM8hfn5+ZQvm2i8u/drsT/qLk4mF0Aklv/jcowDDy+FDMMUNytoUZISQgau6kIgfh+Q9B0grHp4XJ8HjFkIBKwDjPtmuXGPAvW2bdvkvm/ZHONRYrEYeXl5OHbsmMIWkoRoitR8AT6PysRfacV4tPNnkqsFXgoZhvEuZpSkhJCBrugG8N3TQJ1iPgQIq4D4z4AbEcDC4wB/ZK9Xp0ezvltnInt0t6y22Nra4sSJE/Dz8+vuI/s1mvWteRiGQXxWOb6IuouLGWVy51gs4DEvG7wYPAzediZqqiEhpE9VFwIHgtsO0o8ytAJWxyjVslbbrO/IyEgA0g+7yZMnY+nSpViyZIlCOW1tbZiZmcHNzU0uuBOiLhIJg/P/lODzqLtIzquSO8fWZiFszBCsDh6KoZZc9VSQEKIe8fuUC9KAtNzlfcD07b1apR4F6uDgYNm/b926FSEhIXLHCNE0IrEEvyUXYn/0XWSU1sqd47C18Zy/A14IcgbfpG8nixBCNICoQTom3RVJ3wGh7/bqBDOVTSbbunVru+caGxuhpaUFNpuWrhD1EIrE+OlaHr6MyUJBVYPcORMOG0snOGHJBCeYGeqqqYaEELWqKwNSfpSfOKaMhkqgMBlwDOiNWgFQYaCOiYnB+fPnsWHDBvB4PABAeXk5Fi5ciPPnz4PNZmP9+vXYuXOnqh5JSKcEDSJ8dzkXh2KzUV7XJHfO2lgPL0xywXP+DjDU69MFEIQQdWEYoKYIKEqR/6ou6P49G6s7L9MDKvt0+uijj5Ceni43E3zjxo04c+YMhg0bhtraWuzatQs+Pj5YsGCBqh5LSJtKa4Q4FJuD7y/noqaxWe6ck7kB1gQPRZjPEOjpaKuphoSQXscwQFWuYlCuu6/a5+gZd16mB1QWqJOSkjBlyhTZ90KhEBEREZg+fTr++usv1NTUYOTIkfjiiy8oUJNek1dRjwMxdxGRkI+mZoncOQ++MV4KHYrHvPjQpjXQhAwsEjFQfhcovgEUJT8MykKBctcb2QLWnkDORaBZqPxzOaaA7eju1FhpKgvU5eXlGDJkiOz7+Ph4CIVCLFu2DABgZGSExx9/HMePH1fVI8kAJxSJcSNfgNpGEbh6bIy0M4E+u+0W8K3iauyPuovfbxRBLJFfJujvbIaXQoYieLglrYEmZCAQi4D7t+VbycWpgKhOuet5jgB/1IOv0dK10Fwr6bkz70rXSStrzMJez1SmskDN4XBQU1Mj+z4yMhIsFktuFjiXy0VlZaWqHkkGqGKBEF/HZiEiIR+CBpHsuAmHjQW+dlg5yQXWxvoAgOu5Ffg88i4u3FJcTjHV3QovhgzFWEezPqs7IUTFREKgNF0+KJekAeLGzq8FCzAf1ioojwJsvAGDDj4TAtZJk5kos0SLaw2MX6v0q3SXygL1sGHD8Ndff6GxsREsFgvHjh2Dh4eH3NaW9+7dg5WVlaoeSQagtEIBlhy6irLaJoVzggYRvrqYjRNJBXh5siv+SC3C1ewKuTLaWiw8MZKPNSFD4WbTu+NGhBAVa6wFSm7KB+X7twBJc+fXsrQBK3fAZmSroOwlzc/dFcZ8acax9jKTteBaS8v1QRpRlQXqF154AatWrcKwYcOgq6uLnJwc7Nq1S67M9evX4eHhoapHkgGmWCBsN0i3VlbbhK2/pckd09XRwgJfO6yaNBQO5ga9WU1CiCo0VD0YT24VlMsyACiRLFNbVzqe3LqlbOWhui5o/khpxrHLD3J9N7TqCeaYSru7x6/tH7m+W1uxYgUyMjLw9ddfo6GhAS+++CJeeeUV2fn4+HjcuXMHK1euVNUjyQDzdWxWp0H6UVw9HSwc74jlgU6wMtLvpZoRQnqk9v6DceRWQbkyR7lr2QbS7urWQdnSDdDu5bwcxnxpxrHQd6XrpBurpbO7bUf3+e5ZPcr13RVNTU1oaGiAoaEhdHQG55pVyvXdPqFIjHEfXJAbk+6MPlsL0a+HysarCSFqxjDSXNmPLoeqKVTuej0TaWu29SQv82GAVv9fRqm2XN9doaurC11dyvpE2nYjX9ClIA0AQpEE9yrqKVATog4MI20VPxqU68s6vRQAYGD+IBi3aimbOkl3wyFyVB6om5ubcfv2bVRVVUEsFrdZJigoSNWPJf1cbWPXgnSLGmH3riOEdIFEDJRnPhKUbwCNSq5RNh4iP8mLPwowtqWgrCSVBWqGYbBlyxZ8+umncsu02tJeACeDU2OzGLGZSv4V/ggjfcofT4hKiUXSmdYKa5Trlbve1OmR5VCjAK5lr1Z5oFNZoH7//fexY8cO8Hg8LF68GHZ2doN2LJoop7FZjIiEfHwemYkiQRcyAT3A47DhPYT2iSak20RCoDStjTXKykzqZAEWroprlDmmvV7twUZlkfTQoUNwdHREQkICzM3NVXVbMgA1NUvw8/U87Ps7E4XdCNAtFvjZt5upjBDyiMYaoLiNNcqMEj2cWjqApfuDgPygC9vaC9Cj/dr7gsoCdXFxMV588UUK0qRdIrEEv1zPx2d/ZypsNeliaYilE5zw3wsZSi3RsuTqYUWgc29VlZD+raFSOobcOiiXZ0K5Ncp67axRpkmb6qKyQO3s7Izq6t7d6ov0TyKxBCcSC/DfvzOQX/lIgLYwxPoprnhilC20tVgY62jaadITS64ejiz3p9nehABAbemDoJz8MChX5Sp3LduwjTXKI3p/jTLpEpUF6hdffBE7duxAaWkppQklAIBmsQQnkgrw6d+ZuFchPxHFydwA66e44slRttDR1pId97Q1wR/rJ+Hr2GxEXMtDVaslWzwOGwv87LEi0JmCNNEsogagMEnavaxnBNiOUX1SDIaR7pmssEa5SLnr9U3kN6KwGQmYDx0Qa5QHOpUF6jlz5uDixYuYMGECtmzZAh8fHxgbt51r2cHBQVWPJRqoWSzB/5IL8enfGcgplw/QDmYGeHnyMISNGSIXoFuzNtbHO7PcsWHacKQWCFAjFMFIXzpxjMakiUapLgTiH6SZFFY9PK7Pk6aZDFjXvTSTDANUZrexRrlcuesNLKQZtFq3lHmOtByqn1JZZjItLS2wWCwwDNPhVoIsFgvNzUokWB+ABnpmMrGEwW8pBfjvhUxkl8lvN2dvxsHLoa4I8xkCdjsBmpB+pehG5xs3GFpJN27gj2y/jEQszXEtWwp1o4trlO1aZfN68GXEp6CsYTQiM9nixYtpr99BSixhcOpGIfZeyEDWffkAPYTHwcuTh+HpsXYUoMnAUV3YeZAGpOe/e1q6wYMxH2huUlyjXHKzC2uUneUDMn8UYGjR8/chGk1lgfrw4cOquhXpJyQSBqdSi/DfCxnILK2VOzeEx8Ha0GGYN9YOujoUoMkAE79Puf2KgYfBWpst3VdZ6TXKw9tYo8zrSa1JP0UZSUiXSSQM/rxZhL3nM5DxSIDmm+hjbegwzPe1g54OjSeTAUjUIB2T7orStPbPaelI91FuyeLVso+yrmHP6kkGjF4J1HFxcUhOTkZ1dTWMjY0xevRoTJw4sTceRfqQRMLgTFox9pzPwO0S+TSxNsb6WBs6FAv87ClAk4GtMEl+4lhXaOtJg/Cja5R19FRaRTKwqDRQX7p0CcuWLUNmZiYAyE0sc3V1xTfffCMbUCf9B8MwOJNWgj3n7+BWsXyAtjLSw9rQYXiGsoSRgUoilqbVzLsC5F0FsiK7d59Zu4Cxy2iNMukylQXqtLQ0TJ8+HfX19Zg2bRpCQ0PB5/NRXFyMyMhInD17FjNmzMDly5fh4eGhqseSXsQwDM6ll2DP+QykF8kns7E00sOLwUPx/DgHCtBkYGmoAgoSgHtXpMG54DrQVNvpZZ2y9qYgTbpFZYF627ZtaGpqwp9//omZM2fKnXvzzTfx119/4cknn8S2bdtw7NgxVT2W9AKGYXDhn1LsuXAHNwvkA7QFVw9rgl2wcLwjBWjS/zEMUJH1oLX8oMVc+g86T7XJUqJMKxxT6bpmQrpBZYE6KioK8+bNUwjSLWbOnIl58+bhwoULqnokUTGGYRB5uxR7zmfgRr78Gk5zQ12sCR6KheMdwdGlAE36qZYMYi1BOe9K50lEWFrSDSjsxz348geufgnEf6b8c8csVH2mMjJoqCxQCwQCODt3vEmCs7MzBAIlF/GTPsMwDKLu3Mee8xlIyauSO2dmqIvVQS5YFOAIA11aJED6mepC+aBclAJIOkm4pGcC2PsB9uOlQXmIjzQtaGsB64AbEcot0eJaA+PXdv8dyKCnsk9eW1tbXL58ucMyV65cga2traoeSXqIYRjEZJRhz/k7SLpXJXfO1ICNVUFDsTjAEYZ6FKBJPyBuliYPybsK5F2W/lOQ1/l15sMetpTtxwEWIwCtTtb+G/OlGcc6S3rCtZaW604aUUIeUNkn8JNPPolPP/0Umzdvxrvvvgt9/YebJgiFQvz73/9GZGQk1q9fr6pHkm5iGAZxmeX45PwdXM+tlDvHM2DjhUkuWDLBCVwK0EST1VcA+QkPx5cLrnee4UtHHxgy9mFQtvMHDLu5NS9/pDTj2OUHub4bWv2/xDGVdnePX0tBmvSYynJ9l5eXY9y4ccjOzoa5uTn8/f1hbW2NkpISXLt2Dffv34eLiwuuXr0KMzMzVTyy31F3rm+GYRB/Vxqgr+XIB2hjfR28MMkFSyc6wUifZqYSDcMw0nzYrSd9ld3u/DojW8Ch1diytTego6v6+okagMJkoLEa0DOWThyjMWnSikbk+jY3N8fly5fxxhtv4NixY/jzzz9l5/T19bFs2TJ8+OGHgzZIq1tLgL6aXSF33EhfBysDXbAs0AnGFKCJpmiqBwoTgXsPurDzr8q3WNvC0pam2WwJyg7jARO7vqkvmwM4Uo4I0jtU2rdpYWGBQ4cO4cCBA7h165YsM5mbmxvYbAoC6nAlSxqgL2c9EqD1dLA80BnLA51hwqH/NkTNBPnyk76KUzuf9KXPkx9bHuJDaTfJgNQrg5BsNhve3t69cetBTSgS40a+ALWNInD12Bhp1/7+zNdyKvDJuTu4dFd+6QlXTwfLJzphRaALTAwoQBM1EIukWzm2BOW8q0B1QefXWYx4GJTtx0kngXU26YuQAUBlgTo9PR3nz5/Hc889B0tLS4XzpaWlOHbsGKZNmwZ3d3dVPXZQKBYI8XVsFiIS8iFoEMmOm3DYWOBrh5WTXGBtLJ28dz23EnvO38HFjDK5exjqamPpRCe8MMkFPINeGKMjpD115dKu65agXJAINDd0fA3boNWkr/GAnS9gQMNmZHBS2WSyxYsX48KFC8jLy4NWG3/lisViODk5YerUqfjmm29U8ch+pzuTCdIKBVhy6CrKatvfGs+Cq4t3ZrnjZHIhYu7clztnoKuNJROkAdrMkAI06WUSCVB2R74buzyj8+tM7Fu1lv2lCUYo3SYZQDRiMtnFixcxZcqUNoM0AGhra2PKlCmIiYlR1SMHvGKBsNMgDQBltU3YEJEid4zD1sbiCY5YNckF5lzamYf0ksZa6bKolqCcfxUQdpLUSEtHumtUS1C28wdMhvRNfQnpImGzEGnlaagT1cGQbQhPc0/o6+h3fqEKqSxQFxcXw97evsMyQ4YMQVFRkaoeOeB9HZvVaZB+lD5bC4sDnLAqyAUWFKCJKjEMUHWv1djyFWmCEUbS8XUG5vKTvmzH0NIlovFK6kpwNP0oTmSeQHXTwz0PjHWNETYsDIs9F8PKwKpP6qKyQG1oaIjS0o7T6ZWWlsolQiHtE4rEiEjI79I1ejpaOPdqMOzNDHqpVkSjteSxbqyRprzsaUBsbnow6avV2uUaJf7QtnR/ZNLXUODBdreE9Ae3Km5hzbk1KBcq5oGvbqrGkfQjOJV1Cvun7YebmVuv10dlgdrHxwcnT57Erl27wOPxFM5XVlbixIkT8PHx6fYzvvvuO1y8eBHXr19Hamoqmpqa8M0332Dp0qVtlq+urkZ4eDiOHz+O4uJi8Pl8zJ8/H1u3bgWXy1UoL5FIsG/fPnz55ZfIzMwEl8vF1KlTsWPHDri4uHS73t1xI18gN3FMGY3NEhRXCylQDzbVhUD8g+xYwqqHx/V50uxYAeuUy45Ve19x0pe4seNr2IbSiV4tQdlurDQrFyH9VEldSbtBurVyYTnWnFuDiCcier1lrbJAvXbtWoSFhSE0NBR79+5FUFCQ7Fx0dDT+9a9/obKyEuvWrev2MzZt2oTc3FxYWFiAz+cjNze33bJ1dXUIDg5GcnIypk+fjueeew5JSUn46KOPEB0djZiYGIXW/erVq3Hw4EF4enpi/fr1KCwsREREBM6ePYvLly/D1dW123XvqtrGrgXpFjXC7l1H+qmiG+3nmxZWSXd4uhEhzTfNH/nwnEQM3L8lP+mrIqvz5/Ec5buxrTwAbUo1SwaOo+lHOw3SLcqF5TiafhQbfTf2ap1U9n/YnDlz8Oqrr+KTTz5BaGgo9PT0YGNjg+LiYjQ2NoJhGLz++uuYO3dut59x8OBBuLq6wtHRETt37sTbb7/dbtn//Oc/SE5OxptvvomdO3fKjr/11lv48MMP8cknn8hdHxkZiYMHDyIoKAjnzp2Drq50hvTzzz+PWbNmYd26dThz5ky3695VXL3uzXil9J+DSHVh55tCANLz34UBMz8EyjMfTPpKkKa77IgWW5oKs3VgNrJRWfUJ0TTCZiFOZJ7o0jUnMk5g7ei1vTrBTKV/Cu/evRuhoaH4/PPPce3aNeTn54PH42Hy5MlYu3YtHnvssR7df+rUqUqVYxgGBw8eBJfLxebNm+XObd68Gfv27cPBgwflAvVXX30FAHj//fdlQRoAHnvsMYSEhODs2bO4d+8eHBwcevQOyhppZwITDrtL3d88DhveQ0x6sVZEo8TvU26bRQCoKwOOr+i4jKGlfFDmjwbYNKeEDHwiiQgZlRn4M+tPuYljyhA0CZBeng4f6+4P63ZG5X1Wjz/+OB5//HFV37ZLMjIyUFhYiBkzZsDQUD6loKGhISZOnIgzZ84gLy9PNlM9KipKdu5RM2bMQFRUFKKjo7Fo0aI+eQd9tjYW+Nrhq4vZSl+zwM++3UxlZIARNUjHpLuNBVh7yq9dNnWmSV9kwGMYBvdq7iG1LBU3y24itSwVt8pvoUnStRU2rdWKalVYQ0UDcnApI0OaYKG9MWVXV1ecOXMGGRkZsLe3R11dHYqKiuDl5QVtbcVA13Kflvv2lZWTXHAiqUCpJVqWXD2sCHTug1oRtZOIgdSf5SeOKWv0/wHe86RZv/Sp94UMfGUNZUi9nyoLzGnlaV1uNXeGy1acnKxKAzJQCwTShAsmJm1/EBkbG8uV62r5zrRkoHnUzZs34eXlpdQ9AMDaWB9Hlvt3mvTEkquHI8v9ZWlEyQAjbpYuk8qNA3LigNxLQKNyv4sKPOYAQyertn6EaIjaplqkl6fLgvLN8psorivu9Dq2FhvuZu4YYTYCp7JOoaGzFLetmOiawMPcoyfV7tSADNQDiaetCf5YPwlfx2Yj4loeqlqNWfM4bCzws8eKQGcK0gOJuBkoSgFyLkqD873LnU/8UpaesWruQ4iaicQi3Km88zAol91EliALDDrOis0CCy4mLvC08IS3hTe8Lbwx3HQ42A9S1hroGOBI+hGl6xHmGtbrmcoGZKBuaRm31wKurq6WK9fV8p1pL5drey3tzlgb6+OdWe7YMG04UgsEqBGKYKQvnThGY9IDgFgkTVSSEyv9yrsCNHUy5mXlAZTf7Xydc2scU+ksbkL6GQkjQW51rmxMOa0sDf9U/AORpPPJtjaGNvAy94KXhRe8LbzhYe4Brm77XdWLPRfjVNYppZZoWXAssMij9+ctDchA3dmY8qNj2IaGhuDz+cjOzoZYLFYYp+5szLuv6LO14edEOwj1e81NQGGitMWcEyddxyyq6+ACFmDjBThNAhwnAo4TpDtJnXlXuk5aWWMWUupO0i+U1pfKtZTTytJQI6rp9DojXSO5oOxl4QVLA8XdHDtiZWCF/dP2d5r0xIJjgf1T9/dJGtEBG6htbW0RFxeHuro6uZnfdXV1iIuLg7Ozs1xu8uDgYBw7dgxxcXFyyVoAyNZPP3qcEKWIhNKNK3JigdxYIO9ax9s8srQAm5GAU6D0y2F829m+AtZJk5kos0SLaw2MX9v9dyCkl9Q01SCtPE0WlFPLUlFa3/nvtK6WLtzM3WQB2dvCGw5GDmCpYOWCm5kbIp6IkOb6zjgBQdPD3lYTXROEuYZhkcei/pfrW5OwWCysXLkS27Ztw/vvvy+X8OT9999HbW0t3nnnHblrVq1ahWPHjmHz5s1yCU9Onz6NqKgoTJ8+HY6Ojn36HqSfEjUA+dceTPx60GLuqIuapS3tknac+DAwKzMj25gvzTjWWdITrrW0nDJpRAnpRU3iJtyuuC3tvi5PQ2pZKrIFnS9BZYGFobyhci1lV56rbFy5N1gZWGGj70asHb0W6eXpqBXVgsvmwsPco893z1LZftTp6ek4f/48nnvuOVhaKnY1lJaW4tixY5g2bRrc3d279YyDBw8iNjYWAJCamorExERMnDgRw4YNAwAEBgZi5cqVAKQt54kTJyIlJQXTp0+Hj48PEhMTcfbsWfj5+SE6Ohocjnw34AsvvCBLITp79mwUFRXhp59+ApfLRXx8PIYPH96terfoyX6kRIM11UtzZOfESoNzQQIg7mBJnZaOdMMMp0DAMRBwGCfdRKO7qouAyw9yfTdUPjzOMZV2d49fS0Ga9DkJI0GOIEeuC/tW5S00S5o7vZZvyJcLyh7mHjBkG3Z6nSbryee/ygL14sWLceHCBeTl5bW5J7VYLIaTkxOmTp2Kb775plvPWLp0KY4caX823pIlS3D48GHZ9wKBoN1NOYyMFD8YJRIJPvvsszY35Rg6dGi36twaBeoBorFWOuErN04anAsSgY4mtWixpeuWnR60mO38Ab1eWHcpagAKk6UzxPWMpa10GpMmfYBhGJTUl8gCcst6ZWUSgRjrGst1X3taeMKCY9EHte5bGhGonZ2dMWnSJHz77bftllm6dCkuXryIu3fvquKR/Q4F6n6qsUa6RConVhqcC5OAjloF2rqAnd+DruyJ0sCsSzuakYFD0CiQG1e+WXYT9xvud3qdnrYe3M3c5VrL9kb2KhlX1nQ9+fxX2Rh1cXGx3OSstgwZMgRFRUrsZ0uIOgkFDwNzTqx0TTMjbr+8tp40BadToDQ42/lSS5YMGI3iRtyquCUXlHOqczq9ToulhaG8obKA7GXuhWGmw8DWoo2DukplgdrQ0BClpR3P1CstLVXYWpIQtWuoBHLjH3ZlF98AGEn75XU4DwOzU6C0W1tHr+/qS0gvEUvEyBZk42b5wxnYdyruoJnpfFx5CHfIw+5rc094mHvAgE09SaqgskDt4+ODkydPYteuXeDxeArnKysrceLECfj49N4OI2SAETVIu5kba6STrWzHqKalWl8hTcPZslyq+CbQUTYjtoF0JnbLrGxbH0BHt/3yhPQDDMOguK4YN8tvPsyDXZaG+ub6Tq/l6fHkuq89zT1hzjHvg1oPTioL1GvXrkVYWBhCQ0Oxd+9euTXH0dHR+Ne//oXKykqsW7dOVY8kA1V1oXQLx6Tv5Dee0OdJZzEHrOvaLOa6sod5snNigdK0jsvrcqWBuWVWtu1ooBeXgRDSFwSNArnu69SyVKWyb+lr68PD3EPaff3gy45rNyjGlTWFygL1nDlz8Oqrr+KTTz5BaGgo9PT0YGNjg+LiYjQ2NoJhGLz++uuYO3euqh5JBqKiG+2vCxZWSTNx3YiQrgvmj2z7HrWlDyd+5cQB9//p+Jl6xoBDgHTil2MgwB8FaA/IFANkkBA2C2Xjyi2t5Xs19zq9ToulBVeeqywge1t4YyhvKHS06P8HdVLpT3/37t0IDQ3F559/jmvXriE/Px88Hg+TJ0/G2rVr8dhjj6nycWSgqS7sPHkHID3/3dPA6hhpy7qm+OHEr9w4oOxOx9frmzxIxflgVrbNSECLcqaT/kksEeOu4C7SytJkQTmjMkOpcWU7rp1cUHYzc6NxZQ2k8j+THn/8cTz++OOqvi0ZDOL3KZcOE3gYrJuFQEUny/04pg/Hlx0nAtaeFJhJv8QwDIrqimQBObUsFenl6Upty2iqZyo3ruxl4QVT/TZS0xKNQ/0ZRDOIGqRj0l3R3lizgbl8YLbyANpIwkOIqgibhUgrT0OdqA6GbEN4mnuqJM1klbBKbrLXzbKbqBBWdHodR4cDdzN3aVC2lAZnW0NbGlfupyhQE81QmCQ/cawrDC0fBmanQMDSDaAPJNIHSupKpBs3ZJ5AddPDPcONdY0RNiwMiz0XK71xQ0NzA25V3ELq/Yet5fza/E6v02Zpw9XUVa617GLiQuPKA0i3/0tqaWlBS0sL6enpGD58OLS0tJT6a43FYqG5ufOxEzKINDcC97qZrW32J4DvMgrMpM/dqrjV7laI1U3VOJJ+BKeyTmH/tP1wM3OTO98sacbdqrtyk70yqzIh7iixzgP2RvayoOxt4Y0RZiPA0aEEOwNZtwN1UFAQWCwWDAwM5L4npFMMA5SkAVmRQFaUdE2zqPO1m22ycqcgTfpcSV1Jp/sVA0C5sBxrzq3B3tC9cmPL/1T8o9S4spm+mXwebHNP8PR5KnoL0l90O1BHRUV1+D0hcgQF0qDcEpzrOs8L3CmOqXSNMyF97Gj6UaXWIAPSYL3w9MJOy3F0OPA095QLzDaGNtQAIjRGTXqJsFq6XKolOHe0ZIqlLc2PzTDS7SKVNWYh5dQmfU7YLMSJzBM9uocOSweupq5yM7BdTFygTasRSBtUFqhdXFzwyiuvYP369e2W2bdvH3bv3o2srCxVPZZoCrEIKLgO3H3QYs6/1vFGFhbDAZcQwCVUOgFM31i6r/KBIOWWaHGtpfssE9KHmsRN+F/m/+QmjikrwDYAwXbB8DT3hJuZm0pmhZPBQWWBOicnB1VVVR2WqaqqQm5urqoeSdSJYYCyDGlr+W6ktPXcVNN+eUPLB4H5wZeJnWIZY74041hnSU+41tJyXUkjSkg31DbVIuV+Cq6XXEdiaSJult1Eo7ixW/da6L4QQXZBnRck5BF92vUtEAigp0e7DPVbtaUPurIffFUXtF9WhwM4TgCGhkoDs5WncmuZ+SOlGccuP8j13VD58BzHVNrdPX4tBWnSK8oaypBYkojE0kQkliTiduVtSDraSa0LuGyuSu5DBp8eBeqYmBi573NychSOAYBYLEZeXh6+//57DB8+vCePJH2pqV46I7tlAljJzQ4Ks6S7W7mESIOznT/A7mbXnjEfmL4dCH0XKEwGGqul+bhtR9OYNFEZhmGQV5Mnay0nlSYht7rjHj8dLR24m7njduVtNImblH6Wia4JPMw9elplMkj1KFCHhITIZiSyWCwcOXIER44cabMswzBgsVjYuXNnTx5JepNEDBQlPxxnzrsCdPRhZOr0cJzZOQgwMFNtfdgcwDFAtfckg5ZYIsadyjuy1nJiaSLKGso6vMZAxwCjrUZjjNUYjLUeCy8LL3B0OPjo2kc4kt72Z11bwlzDaEyadFuPAvWWLVvAYrHAMAy2bduG4OBghISEKJTT1taGmZkZQkND4e7u3pNHElWryJIG5buRQHZMx9nB9HmAS7A0MLuEAGbOfVNHQrqhUdyI1PupSCpNwvXS60gpTUGtqLbDa8z0zeBj5QMfa+nXCNMRbWb4Wuy5GKeyTim1RMuCY4FFHou6/R6E9ChQh4eHy/49Ojoay5Ytw+LFi3taJ9IeUYM01WZjDaBnJO1q7mpXcH0FkB39sNVc1UFXn7YuYD/uwThzqHT7R1o+QjRUdVM1kkuTZa3lm2U3IZKIOrzG3she1lr2sfKBo7GjUuuWrQyssH/a/k6TnlhwLLB/6n6l04gS0haVTSaLjIxU1a3Io6oLpTtLJX0n3+LV50knVwWsa39ylUgo7cJumZ1dlAKAaf9Z1t7A0BBpi9lhAqBLW94RzVRSVyJtLT8YY86ozADTwe82CyyMMBuBMVZjpC1mK58eBVA3MzdEPBEhzfWdcQKCJoHsnImuCcJcw7DIYxEFadJjLIZhOvjUVl5eXh4yMjIwfvx4WVpRiUSCXbt24bfffgOHw8Grr76K2bNnq+Jx/VJAgHS8NT6+C3mti250vlzJ0Eq6XIk/EpBIpJO+WhKN5MYDHaUqNB4ibS0PfTDOzKUPFaJ5GIZBTnWOrLV8veQ6Cmo7WHUAgK3FhreFtywoj7YaDSNdo16pn7BZiPTydNSKasFlc+Fh7kFj0kROtz7/H1BZi3rz5s34/fffUVxcLDu2Y8cObN26VfZ9dHQ0Ll26BD8/P1U9dmCrLuw8SAPS84dnAU5B0tZzfQcTZHSNAOdJD4Oz+TDKlU00TrOkGbcrbsvNyO5se0cum4vRVqNl3dieFp7Q0+6b5aD6Ovrwsfbpk2eRwUdlgTouLg5Tp04Fm80GIP0L+LPPPoObmxvOnj2L4uJiTJ06Fbt27UJERISqHjuwxe9TLksXIB23vv2H4nEtHcDO7+Hs7CFjAW3KHEs0S0NzA1Lvp+J66XUkliQi5X5Kp5tWWHIsZa3lsdZjMYw3jFJwkgFJZZ/YpaWlcHR0lH2fnJyM+/fvIzw8HHZ2drCzs8PcuXMRHR2tqkcObKIG6Zh0d1iMeJhoxClQOvGMEA0iaBQ8TCxSmoj08nQ0Szre/tbJ2EkWmH2sfGBnZEcbVpBBQWWBWiKRQCJ5mMEnKioKLBYLkydPlh0bMmSIXNc46UBhUsdLpdqz4Cjg8aTKq0NITxTVFuF66XUklSQhsTQRmVWZHZbXYmnBzcxNtlRqjNUYWHAs+qi2hGgWlQVqBwcHXL36cOejkydPgs/nY8SIEbJjxcXF4PF4qnrkwNbYQd7sjuhQilaiXhJGgmxBtmx8ObEkEUV1RR1eo6eth5GWI6VLpazGYpTVKBiyDfuoxoRoNpUF6qeffho7duzAvHnzoK+vj9jYWKxbt06uTHp6OlxcXFT1yIGtu93VesaqrQchnRBJRPin/B8kliTieul1JJcmo6qxqsNrjHWN4WPlgzHWY6QTv8w9wdZm902FCelnVBaoX3vtNZw9exa//vorAGDkyJFyCVFyc3Nx9epVvPXWW6p65MBmO0a6Tror3d8cU2k+bEJ6Ub2oHin3U2St5Rv3b0AoFnZ4jbWBNXysfTDWaix8rH0wlDcUWiwlNmkhhKguUBsbG+Py5cu4eVO6cYO7uzu0teVnYP7666/w9fVV1SMHNjZHmswk/jPlrxmzkDatGMSEzUKklaehTlQHQ7YhPM09VbKWt0JYgaSSJNmM7FsVtyDuaK9xAC4mLnIzsm25tj2uByGDlcrX6Xh5ebV53NHRUW5WOFFCwDrgRoRyS7S41tLtH8mgU1JXIs2OlXkC1U3VsuPGusYIGxaGxZ6Llc6OxTAMCmoL5DauyBZkd3iNDksH7ubuchO/TPVNe/ROhJCHaEGtJjPmSzOOdZb0hGstLUd7NA86typutZtvurqpGkfSj+BU1insn7YfbmZuCmUkjAQZlRnSpCIPWs2l9R3/YcjR4WCk5UhZN7a3hTcM2JRqlpDeotJAXVNTg88++wznz59HYWEhGhsbFcqwWCzcvXtXlY8d2PgjgdUxwOUHub4bKh+e45hKu7vHr6UgPQiV1JV0uikEAJQLy7Hm3BpEPBEBnh4PaeVpstZyUmkSapo6XmFgqmcqlx/bzdwNbC2a+EVIX1FZoL5//z4mTJiAu3fvwtjYGNXV1TAxMUFTUxMaGqQZhmxtbWWZy0gXGPOB6duB0HeBwmSgsVo6u9t2NI1JD2JH048qtc0iIA3Wz//xPKoaq9AoVvwDurUh3CGyGdljrcbC2cSZEosQokYqC9Th4eG4e/cuvv32W/zf//0ftLW18eqrr2LLli24du0aXn75Zejo6ODs2bOqeuTgw+YAjgHqrgXRAMJmIU5knujSNSX1JQrHWGBhmOkwWbYvH2sf2BjaqKqahBAVUFmg/vPPPzFlyhQsXLhQ4Zyfnx9Onz4Nb29vvPfee/jwww9V9VhCBqW08jS5iWPK0mZpw9vCW9ZaHm01GiZ6Jr1QQ0KIqqgsUBcVFWH+/Pmy77W1tWVd3gBgamqKxx57DBERERSoCekmQaMACcUJOJl5slvXfxT8EaY6TlVtpQghvUplgdrExAQikUj2vampKfLz8+XKGBsbo6REsfuNENK2mqYaXC+5jqvFV3Gt+BpuV9wGg+5vIW+mb6bC2hFC+oLKArWLiwtycnJk348ZMwbnzp1DeXk5zM3N0dDQgN9//x0ODg6qeiQhA06dqA6JJYm4VnwNV4uv4p+KfyBhJJ1fqAQTXRN4mHuo5F6EkL6jskA9ffp0fPLJJ6ivr4eBgQFWr16NefPmYdSoUQgICEBiYiJycnKwY8cOVT2SkH6vXlSP5NJkWYs5rTytw6xfXDYXvta+8LXxxa2KWziVdUrpZ4W5hqkkUxkhpG+pLFCvWbMGHh4eskD91FNPYdeuXdi+fTuOHz8ODoeDDRs24PXXX1fVIwnpd4TNQqTcT5EF5tSy1A73YTbQMYCPtQ/8bfzhb+MPNzM3aGtJU/OW1pcivjBeqSVaFhwLLPJYpLL3IIT0HRbDMN0f8FKCWCxGWVkZrKysBv1azIAA6dKq+Ph4NdeE9JUmcRNS7qfIurJv3L8BkUTUbnmODgdjrMbAz8YPfjZ+8DD36DC5SEeZyVpYcCywf+p+jDAb0W4ZQkjv6snnf6+nENXW1oa1tXVvP4YQjSASi3Cz/CauFklbzMn3kztMMKKnrYfRlqPhZ+MHf74/vMy9urTdo5uZGyKeiJDm+s44AUGTQHbORNcEYa5hWOSxSOlc34QQzUO5vgnpgWZJM9LK03Ct+BquFV9DUmkSGpob2i3P1mJjpOVI+Nv4w8/GDyMtR0JPW69HdbAysMJG341YO3ot0svTUSuqBZfNhYe5B41JEzIAdDtQu7i4dOs6yvVN+jOxRIxbFbdwtfgqrhZfRWJJIuqb69str8PSgbelt7TFbOOPUZajei146uvow8fap1fuTQhRn24HaolE0q0x514eEidEpSSMBLcrbstazNdLrqNG1P4mFtosbXiae8oC82ir0bSzFCGkR7odqFuvmSZkoJAwEmRWZUonfxVdRUJJQoepOrVYWnA3c5d1ZftY+8CQbdiHNSaEDHQ0Rk0GNYZhkC3IlnVlJxQnoLKxst3yLLAwwmyErMXsY+0DY13jPqwxIWSw6dNA3dTUBKFQCGNj+mAjnRM2C5FWnoY6UR0M2YbwNPfs8fguwzDIrc6VrWO+Vnyt03XIrqaushazr7UvbWJBCOlTPQrULi4ueOWVV7B+/XrZsTNnzuDMmTP4+OOPFcr/+9//xrZt2yAWt595iZCSuhLpcqPME3Ldzsa6xggbFobFnouVXm7EMAzya/Nl65ivFV1DaUNph9e4mLjIWsy+Nr6UH5sQolY9CtQ5OTmoqqqSO3b58mXs3bu3zUBNSGc6SuBR3VSNI+lHcCrrFPZP2w83M7c271FYWyjXYi6qK+rwmU7GTnKB2YJjoZJ3IYQQVaAxaqIxSupKOs2yBQDlwnKsObcGEU9EwMrACiV1JbLAfLX4KgpqCzq83o5rB3++tCvbz9oP1oaUkIcQorkoUBONcTT9qFJ5q4GHwbpJ0oTc6twOy9oa2soyf/lZ+4HP5auiuoQQ0icoUBONIGwW4kTmiS5dk1GV0eZxKwMr2SYWfjZ+sDOyU0UVCSFELShQE42QVp7W4XrljlhwLGRjzP42/rA3sh/0G8AQQgYOLXVXoLc5OTmBxWK1+RUSEqJQvrGxEdu2bYOrqyv09fVha2uLVatWobS045nCpPvEEjHSytK6de3m8Zvx9/y/8Z+g/2De8HlwMHagIE0IGVB63KL+7rvvcPnyZdn3mZmZAIBZs2YplG0519dMTEzwyiuvKBx3cnKS+14ikWDOnDk4c+YMxo8fj6effhoZGRk4ePAgLly4gMuXL8PS0rJvKj3A5dXkIb4wHpeLLuNK0ZVut6aH8YZRYCaEDGg9DtSZmZltBuC//vqrzfLq+FDl8XgIDw/vtNyRI0dw5swZPPfcc/j+++9ldd2/fz9efPFFbNq0CQcOHOjl2g5MgkYBrhRdQXxRPOIL4zudma0ME10TeJh7qKB2hBCiuXoUqLOzs1VVD43w1VdfAZAmZmn9B8Xq1auxa9cufP/999izZw84HI66qthvNImbkFyaLAvM6eXpYND+hiwjTEdAR0sHaeXKd4GHuYbRNo6EkAGvR4Ha0dFRVfXoVY2NjTh8+DAKCwthbGwMPz8/jBs3Tq6MUCjElStXMGLECIX3YrFYmDZtGg4cOICEhARMmjSpL6vfLzAMgzuVd3C56DLiC+NxveQ6hGJhu+WtDawRYBuAAH4A/Pn+sOBYoLS+FAt+X6DUEi0LjgUWeSxS5SsQQohGGhSzvouLi7Fs2TK5Y35+fvjxxx8xdOhQAMDdu3chkUjg6ura5j1ajmdkZHQaqAMCAto8fvPmTXh5eXW1+hqrpK5E1mK+XHQZFcKKdssasg3hZ+OHAH4AxtuOh7Oxs8IwiJWBFfZP299p0hMLjgX2T92vdBpRQgjpzwZ8oF62bBkmTZoELy8vcLlc3LlzBx9//DGOHj2KKVOmIDU1FUZGRhAIBACkE8/a0rKRSEu5wahOVIdrxddkgTlLkNVuWW2WNkZajpQFZi8LL7C12J0+w83MDRFPREhzfWecgKDp4c/bRNcEYa5hWOSxiII0IWTQGPCBeuvWrXLfjx49Gt9++y0A4OjRo/jqq6+wYcMGlT4zPj6+zePttbQ1VbOkGTfLbsoC8437N9DMNLdb3tnEGeP54xHAD4CfjR+4utxuPdfKwAobfTdi7ei1SC9PR62oFlw2Fx7mHjQmTQgZdAZ8oG7P6tWrcfToUcTFxWHDhg2ylnR7LebqaunyofZa3AMBwzDIqc5BfGE84ovikVCcgFpRbbvlzfTNMI4/DgH8AATYBsDG0Eal9dHX0YePtY9K70kIIf3NoA3UFhbSHZLq6uoASLfs1NLSQkZG22kpW463N4bdX1UIK3C58LJ0ElhRPIrritstq6+tj7HWY6WtZtsAuJq6Qos14HPmEEKIWg3aQH3lyhUAD5OecDgc+Pv74/Lly8jNzZWb+c0wDM6dOwdDQ0P4+vqqo7oApPmw08rTUCeqgyHbEJ7mnl3uChY2C5FYkoj4Iml39q2KW+2WZYEFd3N3WYt5tNVo6Gnr9fQ1CCGEdMGADtS3bt2Cg4MDDAwMFI6/+eabAIDnn39ednzVqlW4fPky3n77bbmEJwcOHEBWVhZWrVqlljXUJXUl0slVmSfkMngZ6xojbFgYFnsubndylYSR4J+Kf2TjzEklSWiSNLX7rCHcIbIW8zibceDp81T9OoQQQrpgQAfqY8eO4eOPP0ZQUBAcHR1haGiIO3fu4M8//4RIJMLbb7+NoKAgWfklS5bgp59+wo8//ojs7GwEBwcjMzMTv/76K5ydnbF9+/Y+f4dbFbfaXa5U3VSNI+lHcCrrFPZP2w83MzcAQEFtgVx6zqrGqnbvb6RrhHE242Rrmu2M7CglJyGEaBAWwzDtp4vq56Kjo/H5558jKSkJJSUlqK+vh4WFBcaNG4eXXnoJ06dPV7imsbERO3fuxNGjR5GXlwczMzM8/vjj2L59O6ytrXtUn5ZZ3+3NCn9USV0Jnjn1jFIJQIzYRgh1CEVyaTLu1dxrt5yOlg5GW46WBWYPcw9oa2kr9wKEEEK6pauf/60N6ECtabr6H+qjax/hSPqRHj93GG+YLDCPtR4LA7ZB5xcRQghRmZ4E6gHd9d2fCZuFOJF5olvXWnIsEWAbgPH88RjPHw9LA9rxixBC+isK1BoqrTytW1s/7pi4A08MfYLGmQkhZICgRbAaqk5U163rePo8CtKEEDKAUKDWUIZsw25dx2V3L20nIYQQzUSBWkN5mnvCWNe4S9eY6JrAw9yjl2pECCFEHShQayh9HX2EDQvr0jVhrmG0aQUhhAwwFKg12GLPxTDXN1eqrAXHAos8FvVyjQghhPQ1CtQazMrACvun7e80WFtwLLB/6n7ao5kQQgYgCtQazs3MDRFPRGCp51KY6MpvsWmia4Klnkvx0+M/YYTZCDXVkBBCSG+iddT9gJWBFTb6bsTa0WuRXp6OWlEtuGwuPMw9aEyaEEIGOArU/Yi+jj58rH3UXQ1CCCF9iLq+CSGEEA1GgZoQQgjRYNT13YeysrJQX18v20WFEELI4HDz5k0YGHRv50IK1H3I1NRU3VXod27evAkA8PLyUnNNBgb6eaoe/UxVa6D+PA0MDLodA2g/aqLRerKHK1FEP0/Vo5+patHPUxGNURNCCCEajAI1IYQQosEoUBNCCCEajAI1IYQQosEoUBNCCCEajGZ9E0IIIRqMWtSEEEKIBqNATQghhGgwCtSEEEKIBqNATQghhGgwCtSEEEKIBqNATQghhGgwCtREI127dg2zZs0Cj8eDoaEhxo8fj4iICHVXq99ycnICi8Vq8yskJETd1dNY3333HVavXg1fX1/o6emBxWLh8OHD7Zavrq7Ghg0b4OjoCD09PTg5OeH1119HbW1t31Vag3Xl5xkeHt7u7yyLxUJOTk6f1l2daJtLonEiIyMxY8YM6Ovr49lnn4WRkRGOHz+OZ555Bnl5edi4caO6q9gvmZiY4JVXXlE47uTk1Od16S82bdqE3NxcWFhYgM/nIzc3t92ydXV1CA4ORnJyMqZPn47nnnsOSUlJ+OijjxAdHY2YmBjo6+v3Ye01T1d+ni2WLFnS5u8oj8dTfQU1FUOIBhGJRMzQoUMZPT09JikpSXa8qqqKGT58OKOrq8vk5OSor4L9lKOjI+Po6KjuavQ7586dk/2+/fvf/2YAMN98802bZbds2cIAYN58802542+++SYDgPnggw96u7oarys/z61btzIAmMjIyL6roIairm+iUf7++2/cvXsXzz//PEaPHi07bmJignfeeQdNTU04cuSI+ipIBpWpU6fC0dGx03IMw+DgwYPgcrnYvHmz3LnNmzeDy+Xi4MGDvVXNfkPZnyeRR13fRKNERUUBAKZPn65wbsaMGQCA6OjovqzSgNHY2IjDhw+jsLAQxsbG8PPzw7hx49RdrQEhIyMDhYWFmDFjBgwNDeXOGRoaYuLEiThz5gzy8vJgb2+vplr2TzExMbhy5Qq0tLTg6uqKqVOngsvlqrtafYoCNdEoGRkZAABXV1eFczY2NuByubIypGuKi4uxbNkyuWN+fn748ccfMXToUDXVamDo6Pe25fiZM2eQkZFBgbqLtm7dKvc9j8fD3r17sXjxYjXVqO9R1zfRKAKBAIC0q7stxsbGsjJEecuWLcOFCxdQUlKCuro6JCUlYdGiRbh27RqmTJmCmpoadVexX1Pm97Z1OdK5UaNG4dChQ8jKykJDQwOys7Px6aefgsViYenSpfjtt9/UXcU+Qy1qQgaBR1slo0ePxrfffgsAOHr0KL766its2LBBHVUjpE1hYWFy3zs5OWHdunVwd3fHtGnTsGnTJjz55JNqql3fohY10SgtLZL2Wh7V1dXttlpI161evRoAEBcXp+aa9G/K/N62Lke6b8qUKRg6dChSU1NlP9eBjgI10SgtY3xtjUMXFxejtra23XFA0nUWFhYApGuASfd19Hvb+jj97qpGy+9tfX29mmvSNyhQE40SHBwMADh79qzCuTNnzsiVIT135coVAJT0pKdcXV1ha2uLuLg4hT966urqEBcXB2dnZ5pIpgJ1dXVIS0uDoaGhLGAPdBSoiUaZMmUKXFxc8MMPPyA5OVl2XCAQ4IMPPoCuru6gmu2pCrdu3Wqz5XHr1i28+eabAIDnn3++r6s1oLBYLKxcuRK1tbV4//335c69//77qK2txQsvvKCm2vU/NTU1uHPnjsLxhoYGvPDCC6ipqcGCBQugozM4plmxGIZh1F0JQlprL4Vobm4uPvroI0oh2kXh4eH4+OOPERQUBEdHRxgaGuLOnTv4888/IRKJ8Pbbb+ODDz5QdzU10sGDBxEbGwsASE1NRWJiIiZOnIhhw4YBAAIDA7Fy5UoA0pbexIkTkZKSgunTp8PHxweJiYk4e/Ys/Pz8EB0dDQ6Ho7Z30QTK/jxzcnLg4uICPz8/uLu7w8bGBiUlJTh//jzy8/Ph7e2NyMhImJubq/N1+o6aM6MR0qYrV64wM2fOZIyNjRkOh8P4+/szx44dU3e1+qWoqChmwYIFjKurK2NsbMzo6OgwNjY2zJw5c5gzZ86ou3oabcmSJQyAdr+WLFkiV76qqop55ZVXGHt7e4bNZjMODg7Mxo0bmerqavW8gIZR9ucpEAiYtWvXMn5+foylpSWjo6PDGBkZMf7+/sx//vMfpr6+Xr0v0seoRU0IIYRoMBqjJoQQQjQYBWpCCCFEg1GgJoQQQjQYBWpCCCFEg1GgJoQQQjQYBWpCCCFEg1GgJoQQQjQYBWpCCCFEg1GgJoQQQjQYBWpCCCFEg1GgJoQQQjQYBWpCSL938uRJsFgsXLp0Sd1V6dTBgwehra2N1NRUdVeF9BMUqAnpBTk5OWCxWHJfurq6sLe3x/PPP48bN26ou4oDhkgkwhtvvIEZM2ZgwoQJbZZJTk7GmjVr4OHhAWNjY+jq6sLGxgbTpk3D7t27cf/+fYVrHv3vp6OjAz6fj7lz5yImJqbb9V2yZAkcHR3x+uuvd/seZHCh3bMI6QU5OTlwdnbG0KFDsXDhQgBAbW0tLl++jLi4OOjp6eHChQuYOHGimmva/x06dAgrVqzAhQsXMHnyZLlzEokEb7zxBnbv3g1tbW0EBQVh5MiRMDQ0RGlpKeLj45GWlgZDQ0Pcvn0bQ4YMkV3LYrFgbm6OdevWAQCEQiGSk5Nx5swZsFgs/PTTT5g/f3636vzpp59i/fr1iI2Npd8B0jn17rJJyMCUnZ3NAGBmzJihcO7dd99lADDBwcF9X7EByNfXl7G3t2ckEonCubfeeosBwPj4+DAZGRltXn/9+nVm6tSpCucBMCNGjFAo/9VXXzEAGCcnp27XubS0lNHR0WEWLlzY7XuQwYMCNSG9oKNAXVxczABgDAwMZMdaAnd+fj6zaNEixtrammGxWExkZKSsTHR0NPP4448z5ubmjK6uLjNs2DDm3XffZerq6hSe8csvvzBBQUGMpaUlo6enx/D5fGbKlCnML7/8Ilfu77//ZmbOnMnw+XxGV1eXsbKyYgIDA5kDBw4ovMuSJUvafNe2/ugIDg5mADANDQ3Mu+++y7i4uDA6OjrM1q1bZWWysrKYFStWMPb29oyuri5jY2PDLFmyhMnJyengJysvNTWVAcC88sorCudu377NaGtrM5aWlkxpaWmn9xKJRArv1VagFovFjKGhIQOAKS0tZaqqqpidO3cyQUFBDJ/PZ9hsNsPn85lFixYxmZmZ7T5v6tSpjL6+PlNTU6PEm5LBTEcdrXhCiLRrtbXy8nIEBATAzMwMzz77LIRCIYyNjQEAX3zxBdauXQsej4cnnngCVlZWSEhIwI4dOxAZGYnIyEjo6urKyr700kvg8/kICwuDubk5iouLcfXqVZw4cQJPP/00AOCPP/7AE088AR6Phzlz5oDP5+P+/ftISUnB0aNHsWrVqh6/49NPP42UlBTMnDkTPB4Pzs7OAIArV65gxowZqKurw+OPPw5XV1fk5OTg+++/x+nTpxEfHw8XF5dO73/hwgUAwPjx4xXOHTlyBGKxGKtXr4alpWWn99LR6frHIYvFwj///IMtW7YgNDQUYWFhMDQ0xK1bt/DDDz/gjz/+QGJiIhwdHRWuDQgIwPnz53Hp0iVMnz69y88mg4i6/1IgZCDqqEW9ZcsWBgATGhoqOwaAAcAsW7aMaW5uliuflpbG6OjoMKNGjWLKysrkzv373/9mADAfffSR7JiPjw+jq6vLlJSUKDy79fVPPfUUA4BJTk7usFxPWtSjR49mysvL5c41NTUxTk5OjJGREZOYmCh37uLFi4y2tjbz+OOPt/msR82fP58B0Ga3dmhoKAOAuXDhglL3ehTaaVEfOnSIAcA4OzszDMMwVVVVCu/IMNLeCi0tLWblypVt3v9///sfA4DZsmVLt+pHBg9qURPSizIzMxEeHg4AqKurw5UrV3Dx4kXo6+tjx44dcmV1dXXxn//8B9ra2nLHDxw4gObmZnz66acwNzeXO/fGG2/g448/xo8//oiNGzfKjrPZbLDZbIX6PHo9AHA4HKXKdcd7770HMzMzuWOnTp1CTk4Otm3bhjFjxsidCwwMxJw5c3Dy5ElUV1fLehTak5+fDwCwtrZWOFdcXAwAsLW1VTgXFRWFqKgouWMhISEICQmRO1ZWVib77ycUCpGSkoK//voLWlpa2LVrFwDAxMSkzbqFhobC09MT58+fb/N8S51b3oGQ9lCgJqQX3b17F++99x4AafC0trbG888/j7feegve3t5yZZ2dnWFhYaFwj8uXLwMAzpw5I+vqbY3NZuPWrVuy75999lm88cYb8PLywvPPP4/Q0FAEBgYqBL1nn30Wv/76K8aPH4/nn38eU6ZMwaRJk9qsQ3f5+/u3+z63b9+WBcHWiouLIZFIcOfOHfj6+nZ4//Lycmhra8PIyKhL9YqKipL9d2nt0UBdXl4uK6etrQ0LCwvMmTMHGzduxKRJk+Tut2fPHly5cgVlZWVobm6WnWsZknhUyx8wZWVlXao7GYTU3aQnZCDqqOu7LQCYoKCgNs8NGzZM1jXe0VcLiUTCfP3114yvry/DYrEYAIyOjg4zZ84cJisrS+7eJ0+eZIKCghhtbW0GAMNisZjJkyczSUlJCu/Sna7vtmZir1y5Uqn3iYqK6vTnNmrUKAYA09TUpHAuJCREqa7vH3/8kQEgN9Gt5b3a6vp+VEREBMNisRgjIyNm3rx5zGuvvcZs2bKF2bp1K+Po6Mi09zGbkpLCAGCeffbZTp9BBjdqUROiIR6dXNaipSVcXV2tVMuRxWJh+fLlWL58OcrLy3Hx4kX8+OOPiIiIQEZGBm7cuCHrXp8zZw7mzJmDmpoaxMXF4ddff8XXX3+NmTNn4tatW+DxeNDSkuZFat1KbCEQCLr8Ti3v8/vvv+Pxxx/v9H060jJJrKKiQqH7e8KECYiKikJkZKTC+mpVCg8Ph76+Pq5fvw5XV1e5c8eOHWv3uoqKCgBQaqIbGdwoMxkhGm7cuHEAHnYZd4W5uTnmzp2Ln376CZMnT0Z6ejoyMzMVyhkZGWHmzJn48ssvsXTpUpSUlODKlSsAAB6PBwAoKChQuC4pKanLdWp5n/j4+C5f+6iW4YPbt28rnFuyZAm0tLTw5Zdf9mr38t27d+Hu7q4QpIuKipCVldXudS11fnQIhJBHUaAmRMO99NJL0NHRwcsvv4x79+4pnK+qqpILmFFRUWAeSTgoEolkLTh9fX0AQExMDMRiscL9SktL5coZGxtjxIgRiI2NlQvyNTU1ePvtt7v8PnPmzIGDgwM+/vjjNlNxikQixMbGKnWv4OBgAJD9UdHa8OHD8cYbb6C0tBSPPfZYm3+gANKfX084OjoiMzMTJSUlsmNCoRAvvvgiRCJRu9e11LnlHQhpD3V9E6LhvLy88Pnnn+PFF1/EiBEjMGvWLAwdOhQ1NTXIyspCdHQ0li5div379wMA5s6dC2NjY4wfPx6Ojo4QiUQ4d+4c0tPTMW/ePNma3vXr16OwsBCBgYFwcnICi8VCbGwsrl69ivHjxyMwMFBWh40bN2LVqlUICAjA/PnzIZFIcPr0afj5+XX5ffT09PDLL7/gscceQ3BwMCZPngxvb2+wWCzk5ubi4sWLMDc3l5sg154pU6bAyMgI586dazN39o4dO9DU1ISPP/4Ybm5uCAoKwqhRo2BgYIDS0lLcuHEDV69eBZfLxejRo7v8LgDw8ssv4+WXX8aYMWMwb948NDc349y5c2AYBqNGjUJKSorCNQzD4MKFC3B3d8fw4cO79VwyiKh5jJyQAak7k8k6Syl69epV5tlnn2VsbW0ZNpvNWFhYMD4+Psxbb73F/PPPP7Jyn3/+OfPkk08yjo6OjL6+PmNubs74+/szX3zxhdykq2PHjjELFixghg4dyhgYGDAmJibMqFGjmA8//LDNbFn79u1jXF1dGTabzTg4ODBbtmxhmpqaOpxM1pH8/HzmX//6F+Pq6sro6ekxxsbGjLu7O7Ny5courX1+8cUXGW1tbaawsLDdMomJicyqVasYNzc3hsvlMmw2m7G2tmYmT57M7Nq1q80151ByMplEImH279/PeHp6Mvr6+oyNjQ2zYsUKprS0tN2fQ1RUFAOA2bNnj9LvSQYv2pSDENKv3b59G15eXggPD8e7776r7uooZeHChTh9+jTu3r0rmwNASHtojJoQ0q+NGDECK1euxCeffIKamhp1V6dTd+7cwbFjx7Bp0yYK0kQpNEZNCOn33nvvPVhbWyMnJ0fjZ1Hn5+dj69atWLt2rbqrQvoJ6vomhBBCNBh1fRNCCCEajAI1IYQQosEoUBNCCCEajAI1IYQQosEoUBNCCCEajAI1IYQQosEoUBNCCCEajAI1IYQQosEoUBNCCCEajAI1IYQQosEoUBNCCCEajAI1IYQQosEoUBNCCCEajAI1IYQQosH+H+/pKz6nRczUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x364 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.6, 2.6), dpi=140)\n",
    "\n",
    "for col in df:\n",
    "    if not col.startswith('c'):\n",
    "        continue\n",
    "    ax.plot(df.pressure, df[col], 'o-', label=col)\n",
    "\n",
    "ax.set_xlabel('Pressure (GPa)')\n",
    "ax.set_ylabel('Elastic constant (GPa)')\n",
    "ax.legend(frameon=False)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcb3b571-11ed-4605-a383-bc80c71371e5",
   "metadata": {},
   "source": [
    "## CsPbI<sub>3</sub>\n",
    "\n",
    "In materials with internal degrees of freedom one can distinguish the so-called relaxed $c_{ij}$ and clamped elastic constants $c_{ij}^0$.\n",
    "In the former case, the ionic positions are allowed to relax after application of a macroscopic strain, whereas in the latter the relative coordinates are kept fixed.\n",
    "\n",
    "To illustrate this effect, we here consider the elastic stiffness tensors in orthorhombic structure of CsPbI<sub>3</sub>, which we describe using a neuroevolution potential (NEP) taken from [Fransson *et al.*](https://pubs.acs.org/doi/10.1021/acs.jpcc.3c01542)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "37012612-3abc-4d5c-99db-489c8d57edd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "structure = read('CsPbI3-orthorhombic-Pnma.xyz')\n",
    "calculator = CPUNEP('nep-CsPbI3-SCAN.txt')\n",
    "structure.calc = calculator\n",
    "relax_structure(structure, fmax=0.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e961bfe6-5d46-40e0-874b-099d51d6fa03",
   "metadata": {},
   "source": [
    "Next we can compute the relaxed (`clamped=True`; default) and clamped (`clamped=False`) elastic stiffness tensors.\n",
    "Since this is a rather soft material, we tighten the condition on the convergence of the forces (`fmax=1e-5`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f8dc5654-9c61-4c1c-8c81-7260033d6d1f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relaxed elastic constants\n",
      "[[16.7 14.9 10.1 -0.  -0.   0. ]\n",
      " [14.9 25.3  9.5  0.  -0.   0. ]\n",
      " [10.1  9.5 31.1 -0.   0.  -0. ]\n",
      " [-0.   0.  -0.   6.1  0.  -0. ]\n",
      " [-0.  -0.   0.   0.   4.6 -0. ]\n",
      " [ 0.   0.  -0.  -0.  -0.  11.4]]\n",
      "\n",
      "Clamped elastic constants\n",
      "[[31.1 17.7  5.9  0.  -0.  -0. ]\n",
      " [17.7 34.9  5.2  0.  -0.   0. ]\n",
      " [ 5.9  5.2 46.8 -0.   0.  -0. ]\n",
      " [ 0.   0.  -0.   7.  -0.  -0. ]\n",
      " [-0.  -0.   0.  -0.   6.5  0. ]\n",
      " [-0.   0.  -0.  -0.   0.  18.4]]\n"
     ]
    }
   ],
   "source": [
    "cij_rlx = get_elastic_stiffness_tensor(structure, fmax=1e-5)\n",
    "cij_clamped = get_elastic_stiffness_tensor(structure, clamped=True)\n",
    "\n",
    "with np.printoptions(precision=1, suppress=True):\n",
    "    print('Relaxed elastic constants')\n",
    "    print(cij_rlx)\n",
    "    print('')\n",
    "    print('Clamped elastic constants')\n",
    "    print(cij_clamped)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e684b03-baf3-4589-9799-7636ed7ad7de",
   "metadata": {},
   "source": [
    "Evidently in this material there is a pronounced difference between the relaxed and clamped elastic constants, which is related to the rather soft interactions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
